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PLASMA  THEORY  and  SIMULATION 

A.  LOWER  HYBRID  DRIFT  INSTABILITY  SIMULATIONS  USING  ESI  HYBRID  CODE 

Yu-Jiuan  Chen  (Or.  Bruce  Cohen  (LLL)  and  Prof.  C.  K.  Birdsall) 


In  Che  last  QPR,  it  was  shown  chat  our  simulations  exhibited 

strong  ion  trapping  at  saturation  in  the  lower-drift-velocity  regime  with 

v_ < v  .  as  our  model  had  assumed  a  constant  density  and  v_  drift.  The 
E  ti  E 

typical  total  ion  momentum  (Fig.  la)  and  ion  kinetic  energy  (Fig.  lb)  are 
exponentially  growing  in  time  in  the  linear  regime.  In  Fig.  la,  the 
net  ion  momentum  increases  in  the  negative  x  direction  which  is  the  wave 
propagation  direction.  Recall  that  ions  in  our  model  are  initially  in  an 
electrostatic,  steady  equilibrium  state  with  the  ion  diamagnetic  drift 
velocity  cancelled  by  the  g*8  drift  velocity,  i.e., 


v 


E 


*i 


Vci 


(1) 


or  Vj.  *Vg-vA.*0.,  where  v^.  is  the  ion  drift  velocity.  Also,  ions  are 
shielded  by  the  strongly  magnetized  electrons  which  are  represented  by  a 
linear  electron  susceptibility,  x  (i(3/3t)-kv  ,  v  (T  )),  i.e., 


xe(i  It  '  kvde’v*e(Te))k  \ 


^P, 


(2) 


where  v 
vde  IS 


d  2 

*T/ma)  \  Zn  n|  is  the  electron  diamagnetic  drift  and  T  -m  v 
*e  e  e  ce  1  dy  '  3  e  e  te 

the  electron  drift  velocity  given  by 
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From  the  momentum  conservation  law, 

nm  v  (0)  ■  P.(t)  +  nm  v  (:)  -  P.  (0)  and  P.  (0)  *  0  ,  (4) 

e  de  i  e  de  i  i 

the  increase  of  ion  momentum  in  the  direction  of  wave  propagation  implies 

that  current  relaxation  (v  ,{t)  -*-0)  occurs  during  the  growth  of  the  instab- 

d 

ility,  where  v  .  =  v  ,  -v..  is  the  cross-field  current  velocity  or  the  elec- 
1  d  de  d i 

tron  velocity  in  the  ion  frame,  v^  is  no  longer  constant  in  time  in  the 
late  stage  of  growth.  Similarly,  it  is  noted  that  the  electron  random 
kinetic  energy  must  adjust  in  response  to  the  lower-hybrid  drift  instability 
according  to  the  energy  conservation  law, 

W.  (0)  +  nTe(0)/2  +  nmev^(0)/2  +£(0) 

*  W.  (t)  +  nT  (t)/2  +nm  v*  (t)/2  +£(t)  (5) 

i  e  e  de 

and  the  change  of  ion  kinetic  energy  W.  (Fig.  lb)  and  the  electric  field 
energy  ^(t).  Therefore,  we  model  the  implied  changes  in  the  electron 
kinetic  energy  and  momentum  by  allowing  the  electron  temperature  and 
drift  velocity  vde  to  adjust  in  the  linear  electron  susceptibility. 

To  incorporate  the  effect  of  anomalous  transport  in  our  simula¬ 
tions,  the  1  inear  electron  susceptibi 1 i ty  is  altered  to  accomodate  vde  and 

T  slowly  changing  in  time.  If  time  nAt<nthAt,  both  v,  and  T  are  con- 
e  06  6 

stant  in  time.  Therefore,  using  a  large  value  for  nth  allows  us  to  simu¬ 
late  the  linear  regime  of  the  instability  without  considering  anomalous 
transport. 

If  time  tanAt>nthAt,  we  use  the  predictor-corrector  scheme. 

Set  predicted  quantities  and  T^  as 
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dep 


n-1 

/de 


and 


rn-1 


ep 


(6) 


8y  solving  Eq.  (2),  we  get  the  predicted  electric  field  En  and  field 

P 

energy.  Push  the  particle  velocities  from  vn  ^  to  Vp^  by  using  e” 
and  calculate  the  predicted  total  ion  kinetic  energy  W1^  and  ion  mom¬ 
entum  P?  by 
'P  7 


W. 


iP 


N  , 

».  r  vn-*vn+* 
•  >1  j  jp 


(7) 


and 


Pn  =  (pn-i  +  Pn+i)/2 
ip  \  i  ip  ) 


(8) 


N 


where  P?p  *  m.  21  vjp  *  Now  we  use  the  corrector  portion  of  the 
algorithm.  Eq.  (4)  yields 


de 


v°  +  (P?  -P?)/nm 
de  ip  i  e 


We  substitute  v^  into  Eq.  (5)  and  obtain 


(9) 


nT  nT  nm  v , 

-2.  -  -2-  +  w?  +£°  +  -  w? 

2  2  1  2  'P  P 


n 

nm  v, 
e  de 


(10) 


Substituting  v^e  and  into  Eq.  (2),  we  can  solve  the  electric  field  En 
and  the  field  energy £n. 

Comparison  of  simulation  results  for  time  varying  v^  and  Tg 
with  that  for  constant  v^ft  and  is  given  in  Figs.  2  and  3.  Parameters 
used  are  m./m  ■  1 600 .  ,  /cu^  »  1 . ,  T  *0.,  v  .  *  (T./m. ) ^  ■  0. 141421 ,  and 

I  6  p©  CC  C  1 1  II 

k\p  *  0.707.  Fig.  2  for  Vg"0.12  and  nthAt*100  shows  that  current  relax 
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ation  stabilizes  the  lower-hybrid  drift  instability  (Fig.  2a);  otherwise 
ion  trapping  stabilizes  the  instability  (Fig.  2b  and  Ref.  1).  By  assum¬ 
ing  that  the  system  stabilizes  via  current  relaxation  ,  the  saturation 
level  of  the  fluctuation  spectrum  was  estimated  at 


by  Davidson  (Ref  2),  and  Oavidson  and  Gladd  (Ref.  3).  The  corresponding 
change  in  the  ion  temperature  after  saturation  was  given  as 


AT. 

— -  2 

Ti  ~ 


03) 


as  consequence  of  Eqs.  (12),  (13),  momentum  and  energy  conservation,  Eqs. 

(4)  and  (5),  and  relaxing  the  relative  drift  v.  to  zero.  For  the  para- 

d 

meters  we  used,  AT  /T.  and  AT./T.  are  about  10  The  simulated  satura- 
e  i  it 

tion  level s  for  both  constant  and  time  varying  v^e  and  Tg  compared  with 
theories  (Refs.  1,  2  and  3)  are  given  in  Fig.  3,  indicating  fairly  good 
agreement.  For  cases  in  which  v_ <  v  .  and  v,  is  allowed  to  relax  in 
time,  the  saturation  level  of  £s/nT.  is  so  small  that  it  is  masked  by 
the  thermal  f  1  ucatuat ions ,  whose  level  is  equal  to  MO  in  our  simula¬ 
tions.  Thus,  no  data  appears  on  the  curve  labeled  (b)  in  Fig.  3  for 

v./v..  <  0.8. 

E  1 1 


Comparison  of  saturation  field  energies  as  functions  of  v^/v^... 

(a)  For  constant  v_,  and  T  ,  the  lower-hybrid  drift  instability 
de  e 

is  saturated  by  ion  trapping  [1].  (b)  For  time  varying  vdg  and 

T  ,  saturation  of  the  instability  is  achieved  via  current  relaxa 
tion.  The  simulation  data  (dots  in  (a)  and  crosses  in  (b))  are 
in  agreement  with  the  theory  shown  by  lines. 
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We  would  like  to  thank  Or.  W.  Nevins  for  several  helpful  discussions. 
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B.  CONTROL  OF  UNWANTED  BEAMING  INSTABILITY 
(Prof.  C.  K.  Birdsall) 

Next  quarter  it  is  planned  to  look  into  instabilities  which  occur 
in  magnetized  plasmas  with  a  Maxwellian  or  other  f(v  )  made  up  of  rings  in 
v  space  at  t=0.  The.  ring  may  be  real  (multivelocity  neutral  injection) 
or  the  result  of  loading  for  simulation.  The  multi-ring  dispersion  rela¬ 
tion  will  be  solved  for  complex  and  real  k.  Instabilities  are  expected 

even  with  a  Maxweilian  envelope.  Simulations  may  be  done  to  find  sa¬ 
turation  levels  and  detailed  ring-ring  interaction  in  order  to  aid  in  find¬ 
ing  means  to  control  these  physical  (but  usually  unwanted)  instabilities. 

C.  FIELD  REVERSED  PLASMA  SIMULATIONS,  QUAS (NEUTRAL,  IN  2d 
Douglas  Harned  (Prof.  C.  K.  Birdsall,  Alex  Friedman) 

Progress  has  been  made  toward  solving  several  problems  that  have 
arisen  in  our  two-dimensional  quasineutral  code  AQUARIUS.  Much  of  the  pro¬ 
gress  has  consisted  of  debugging  and  adding  improved  diagnostics.  Some  of 
the  difficulties  originally  encountered  (e.g.,  noise,  growing  E,J)  in  set¬ 
ting  up  cylindrical  problems  in  cartesian  coordinates  have  been  reduced  by 


11 
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making  the  initialization  more  self-consistent.  We  are  currently  working 
on  boundary  conditions  for  a  conducting  wall  boundary  which  are  consistent 
with  the  quasineutral  equations. 

D.  DISTRIBUTION  FUNCTION  AND  CHARGE  AND  CURRENT  DENSITIES  IN 

LINEARIZED  SIMULATIONS  WITH  CYLINDRICAL  COORDINATES 

Alex  Friedman 

In  the  following  we  derive  appropriate  expressions  for  the  zero- 

and  first-order  distribution  function  associated  with  the  kth  simulation 

superparticle.  From  these  expressions  we  calculate  charge  and  current 

densities.  We  then  compare  these  results  with  the  work  of  Bobroff 1 ,  and 

2 

make  application  to  the  RINGHYBRID  3d  code  ;  axial  displacements  and 
current  densities  are  omitted  here  as  they  are  straightforward.  The 
development  presented  here  may  be  contrasted  with  that  of  Cohen  et  al., 
elsewhere  in  this  QPR. 

The  distribution  function  f.  (r,9,v  ,v  )  can  be  written 

k  r  0 

i  r2ff  5  [e-0  (d>)  ] 

fk  “  17  JQ  d*  str“r|t(*)]  - - -  5  [Y"¥k(r-Ere  *.*)].  (0 

where  4  parameterizes  the  unperturbed  angular  location  of  the  element  of 
the  ring  between  $  and  (j>+d$;  and  9  are  both  defined  relative  to  the  labor¬ 
atory  x-axis,  and  9^($)  is  the  angular  location  (in  9)  of  the  point  whose 
unperturbed  location  is  $.  S  is  the  finite-size-particle  analogue  of  the 
<5-function.  The  particle  velocity,  v^,  is  a  vector  field  since  it  differs 
for  different  points  within  the  superparticle.  In  particular,  to  ensure  a 
rigidly  rotating  zero-order  superparticle  (no  shear  to  the  motion),  we  must 


have 
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O  ,  0\  r  _  • 

v9k(rk)  “  =  r9 


(2) 


where  r°  is  the  radius  of  the  "center"  of  the  superparticle.  Note  that  the 
k 

flow  velocity  as  a  point  in  space  (r,9)  arises  from  the  motion  of  an  ele¬ 
ment  whose  unperturbed  location  is  (r-e^e 1  ^  ,  9-e^e 1  ^/r) ,  where  e^e'2,1^ 
is  the  perturbation  from  circularity.  For  rigid  rotation 

£9k  3  C9k^rk^r/,rk^  ’  and  we  have  t*le  r®lati°n: 


eek^rk'  U* 

+  - —  e 


(3) 


Other  expansions  needed  are: 


rk(lf,)  *  rk  +  £  rke 

vrk(*}  “  vrk  +  vrk(^  (4) 


v9k(r^} 


o  i  o, 

vV 


r  .  1  ,  o  .  r 

T  +  V9k(rk’^  ~ 
rk  rk 


where 


VU)  *  ("rk  '  e9kl)e'^ 

(5) 

**<♦>  -  ^k^rk^*  • 


(The  arguments  (r) ,  (r°)  have  been  omitted  in  the  above  for  the  sake  of  sim¬ 
plifying  the  notation;  in  the  final  expressions  their  locations  will  be  evi¬ 
dent.  ) 
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The  zero-order  distribution  function  is: 


rjr  S(r-r°)  6[v-v°(r,9)]  . 


Integration  over  velocities  yields  the  zero-order  density: 


o  1  c  /  o, 

n,  »  » —  S (r-r  )  , 

k  2wr  k'  ’ 


which  satisfies  the  normalization  condition  jdr  2rrr  n^(r)  ■  1  if 

dr  S(r-r°)*1,  as  for  a  6-function  (the  limits  of  integration  are 

0  and  r  ,  the  radius  of  the  outer  system  boundary).  For  simplicity, 
max 

so  that  we  may  freely  interchange  3S/3r°  and  -3S/3r,  we  demand  S(r-r°) 
be  a  symmetric  function  of  (r-r°) .  The  zero-order  azimuthal  current 
density  is: 


2^S(r'rk)  V9 


<(r)  v°k(r) 


so  that  the  algorithm  used  to  accumulate  jQk  on  a  grid  di ffers  from  that 
used  for  n£  if  there  is  no  shear;  for  the  former  one  divides  by  the  parti¬ 
cle  radius  r°,  while  for  the  latter  division  is  by  the  cel  1  radius  r.  Most 
codes  with  which  the  author  is  familiar  use  the  same  algorithm  for  weighting 
n°  and  j°k>  and  thus  imply  shear  within  the  motion  of  each  superparticle; 
this  may  have  particular  significance  in  gridded  ("spokes  and  wheels")  r-6 
simulations  where  shear  implies  deformation  of  the  particle,  which  is  not 
permitted  by  the  9-weighting  algorithm. 


Tt 


-  14  - 


The  first-order  distribution  function  may  be  written  as 


k  *  2^7  /  S[r-rk"erkeU^  5f9-^s9ke'^/rk] 

'A 


5K(^rke'4  -  9-eflkea*/r°)l  -  f°  . 


9k 


(9) 


Integration  over  velocities  yields  the  first-order  density: 

i  40 


1  1  f  3  » /  0\  i  J 

”  '  2^(‘£r  5TS(r'r  )e 


J  d$  5 (0-$-e9e' £9/r°)  -  !  S(r-P°)l 


do) 


where  subscripts  k  have  been  omitted  for  simplicity.  Expansion  of  the 


term  in  square  brackets  yields 


dd>  (ea/r0)e' ^  96  (9~d> )/ 3<}> ,  and  Integra- 


I  0  q 

tion  by  parts  yields  -  i4e  e  /r  .  Use  of  Eq.  (7)  yields: 


_r  _3_ 
r  3r 


i  4e„ 


r  - 


o  i  49 
n  e 


(11) 


which  corresponds  to  Bobroff's  expression  [17]  erian t)  »  _y. (p^r^ ) t 
|  £0 

where  r 1  is  our  se  and  p  our  n.  With  this  information  the  first-order 
distribution  function  can  be  simplified.  We  have 


<  _  l  „  I  *,«  ,  Ox 

vk(r-sre  Y  ,  9-e0e  Vr  ) 


i  4$ 

o,  e9® 

4(r'9)  — — 


iv’fp.S)  -  ■ 

(12) 


and  using  the  relations  3r/39*9,  39/30  * -r  (carets  denote  unit  vectors), 
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vj(r,8)  -  -LT-  {v°§(8)  -  v°  ( r ) f ( 8 ) ] 


-  €  e'^ee  +  v^(r,9) 
r  -k 


y°(r,9)  -  (ee1 i$-7)y°(r,9)  +v'(r,9)  , 


so  that 


5  [v  -  v  ]  -  5  [v  -  v°] 


[(MU*.7)vJ  -  y’]-i-4[y-yj]  , 


and  we  obtain: 


f1  *  n'tS  [y-y°(r,9)  ]  +  n°  [  (ee '  •  7)y° (r , 8)  -  y1  (r  ,9)] &  [v-y°(r,9)] 


_Or/ _ iiLS  „,  _o. 


Note  that  Eq.  (13)  might  have  been  written  down  from  first  principles. 
From  this  the  first-order  current  may  be  calculated: 

i 1  *  /d3y  -  f1 

-  n1y°(r,9)  +  n°  y  [  (ce' 110  •V)v°-y1]  •  5  [y-y°]  ,  (16) 

and  integration  by  parts  yields: 


.1  1  o 


j'  ■  n'yu(r,9)  +  nu[v'(r,9)  -  (ee ' * V) v° ( r , 9 )] 


This  expression  is  in  agreement  with  8ogroff's  Eq.  [B63] , 


{\  *  °°C3T-1  *  (V7)-r1  ‘  (V7)*o]”  V^o-I*  *  [B63] 


if  we  identify  3r,/3t+(v  *V)r  ,-v,  using  Bobroff's  [ 1 0 1  and  our  (11)-  Me 

-  i  “O  “I  “I 


that,  explicitly, 


dee  .  •  ^  i 

=  [-v°(e0/r°)0  +  (er-iZ0er)r  +  (eQ- i  Jt6e0) 9]e 


(vo-v)£eli9  =  [(iJ,9er-0£9)r  +  ( i  J,9£9+9er+v°e0/rO)  9}e‘ 


,  i  Jt8  _i  o 
(-£e  *7)v 


r  o,o?  •  ^  7  I  Jt0 

[e00r-e0vr/r  0-Er00]e 


We  may  express  the  first-order  current  in  component  form  by: 


-1  -  -1-  e-iUa4s(r-rV£v°is(r-r0)  eil6 


r  2irr 


r  9  ro 


.1  1  J /•  r  .  •  \  . /  o>  o  3  o> l  i  £9 

J9  "  ‘  'Z0e9  S(r'r  5  ’  Ve  ¥7  S(r’r  YC 


Note  the  presence  of  r  and  r°  in  the  demominators  of  j  ,  j!,  respectively 

r  y 

(An  alternative  derivation  of  uses  the  identity: 


r2ir  (£  e  "T  In 

f0  •sta-*]  |«ev-v >] 


/.2ir  ,  f£fleiJ,l<l  •) 

-fo  d*  aV  { — ^ —  5[-'-  ]j 


^XU05[v-v°]  -  4-«[x-v°(9)1 
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The  first  term  gives  part  of  n^v°;  the  second  gives  part  of  (-ee ' v°, 
s  i  nee 


_3_ 

30 


<5  [v-v°] 


_3_ 

3v 


6 [v-y°] 


(23) 


through  use  of  the  "Chain  Rule"  of  calculus. 

Using  the  results  of  Bobroff  it  is  not  difficult  to  show  that 
charge  continuity,  7*]1  + 3p 1 /3 1  -  0,  is  safisfied.  An  alternative  expres¬ 
sion  for  the  first-order  current  is  Bobroff's  [857]  : 

j1  *  3 (Por 1 )/3t  +  7  X  (Pq^1  X v°)  .  [B57] 


Taking  the  divergence  of  this  expression  gives  7- j 1  -  3 (7*pQr 1 ) /3t  which 
i  S  -3p  Vst  by  (11). 

As  a  check  on  the  calculation,  the  continuity  equation  can  be 
vefified  directly  from  our  final  expressions  for  n1,  j'.  Since  e  is  de¬ 
fined  as  a  moving  point  which  at  the  old  time  level  is  at  angle  a*0  and 

•  1  new 

at  the  new  time  level  is  at  a*6At,  we  must  rotate  n  back  to  the 
x-axis  through  an  angle  £0At: 
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Taking  the  divergence  of  expressions  (21)  yields: 
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(5)  {3b}  {6}  {4} 

so  that  the  two  expressions  are  in  agreement.  An  alternate  derivation 
of  3nV3t  involves  writing  it  as  -7*  (3n°/3t) r1  -  V*n°3r1/3t  and  recogniz¬ 
ing  that  3n°/3t  *  (-v°/2irr)  3S/3r . 

Another  alternate  derivation  uses  the  relation 
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The  vr3/3r  and  v^/r  terms  combine  to  yield  a  (v^/r) (3/3r) r  term,  and  there 
are  factors  of  1/2irr  in  n1.  Also,  n1  decreases  (even  at  the  center  of  the 
ring)  as  the  ring  moves  outward  (v°>0)  because  of  the  division  by  the 
increasing  value  of  r(ring  center),  hence  the  v°/r  term  in  dnVdt. 
**********  **********  **********  ********** 

As  an  example,  we  consider  the  displacement  off-axis  (£■!)  of  an 


infinite  current  layer  in  a  uniform  field.  The  equations  of  motion  are: 


-  19  - 


e  -  e  92  -  2e  0  *  -  9  (e  +  e  0) 

r  r  9  9  r 

e0  -  ^9  +  2er9  *  9(er  '  £99) 


or 


9e. 


-  9e 


Furthermore,  if  the  vectorial  first  order  displacement  of  a  point  is  to  re- 

•  ,  *  *  .  * 

main  constant  as  the  point  gyrates,  we  must  have  er“'9er>  e0»i9eQ. 
Compatibility  of  these  with  the  above  requires 

z  *  ie  ,  e  »  i9e  ,  d/dt  M  i9 

9  r  r  r 

where  d/dt  is  the  rate  of  change  of  a  quantity  (such  as  zf)  along  a  particle 
orbi t. 


If  motion  «e'iwt  is  assumed,  u  -  -9,  i ndi cati ng  that  the  comoving  coordinate  axes 
rotate  in  the  +9  direction,  so  components  of  a  constant  vector  in  comoving 
coordinates  rotate  in  the  -9  direction. 
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ABSTRACT 

Conservation  laws  are  derived  which  are  appropriate  for  use  in 
linearized  particle  simulation  codes.  A  continuity  equation  describing 
mass  or  charge  conservation  and  an  extended  Poynting  theorem  describing 
energy  flow  are  obtained  for  general  curvilinear  coordinate  systems. 

These  conservation  equations  are  expressed  in  terms  of  quantities  readily 
calculable  in  a  particle  code,  serve  to  quantitatively  assess  the  validity 
and  precision  of  code  results,  and  provide  valuable  insight  into  the 
physics  being  studied.  A  conserved  quantity  related  to  the  system  energy, 
rather  than  a  Poynting  theorem,  is  derived  from  a  Lagrangian  formulation. 
It  Is  bilinear  in  first  order  quantities  but  is  not  attractive  for 
Implementation  as  a  diagnostic  in  a  linearized  particle  simulation. 
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I.  INTRODUCTION 

Particle  simulation  provides  a  powerful  tool  for  assessing  the 
stability  of  complex  plasma  equilibria,  wherein  kinetic  effects  play  an 
essential  role.  Of  special  interest  is  the  recent  use  of  linearized 
particle  codes1"^  to  study  the  stability  of  equilibria  for  which 
analytical  methods  are  not  available.  These  codes  solve  the  linearized 
equations  of  motion  for  small  deviations  of  the  particle  orbits  and 
electromagnetic  fields  from  their  equilibrium  values.  Compared  to  fully 
nonlinear  simulations,  linearized  simulations  possess  certain  advantages. 

By  means  of  linearization  and  Fourier  analysis  the  computer  economics  of 
two-  and  three-dimensional  simulation  are  greatly  improved. 

Furthermore,  highly  accurate  determination  of  linear  growth  rates  is 
possible;  this  is  often  very  difficult  in  fully  nonlinear  simulations 
because  of  high  thermal  fluctuation  levels  and  relatively  short  growth 
periods  terminated  by  nonlinear  saturation. 

Linearized  codes  are  particularly  relevant  to  situations  in  which  the 
dimensionless  ratio  of  a  characteristic  parameter  associated  with  particle 
dynamics  to  a  characteristic  parameter  associated  with  either  the 
equilibrium  configuration  or  a  coherent  oscillation  of  the  plasma  is  of 
order  unity,  so  that  an  analytical  perturbation  expansion  is  not  possible. 
Examples  of  importance  in  the  study  of  fusion  plasmas  occur  when  the 
characteristic  dimension  of  a  particle  trajectory,  for  example,  the  electron 
or  ion  larmor  radius  in  a  magnetized  plasma,  is  comparable  to  either  the 
length  over  which  the  plasma  equilibrium  varies  or  the  wavelength  of  an 
unstable  mode. 
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In  this  article  we  investigate  certain  conservation  properties  of 
linearized  particle  codes.  Conservation  laws  are  valuable  in  at  least  two 
different  ways:  they  provide  quantitative  means  of  monitoring  the  accuracy 
and  validity  of  a  simulation  and  can  furnish  considerable  insight  into  the 
details  of,  for  example,  particle  and  energy  transport  accompanying  the 
linear  phase  of  an  instability.  These  conservation  laws  act  as  the  basis 
for  code  diagnostics  that  only  require  simple  computations  involving 
readily  available  simulation  quantities.  Particular  attention  will  be 
given  to  the  form  these  expressions  take  when  curvilinear  coordinates  are 
used. 

It  is  not  a  trivial  task  to  develop  conservation  laws  for  linearized 

particle  or  fluid  models.  This  is  because,  by  definition,  in  a  linearized 

system  there  are  only  quantities  that  are  linear  in  a  perturbation  and  no 

terms  of  second  or  higher  order.  The  difficulty  this  poses  can  be  seen  by 

considering  as  an  example  particle  kinetic  energy.  With  particle  velocity 
-*■ 

given  by  v  ®  vQ  +  ev^  +  e  Vg  +  ...»  where  e  is  a  small  parameter 
measuring  deviation  from  equilibrium,  the  kinetic  energy  is  K  =  (m/2) 

*[vg  +  2evg  •  v-^  +  (v^  +  2  Vg  •  Vg)  +  ...].  Typically,  the  first 

order  terms  (proportional  to  e)  will  vanish  when  K  is  summed  over 
particles.  The  second  order  kinetic  energy,  when  summed  c/er  all 
particles,  plus  the  analogous  electric  and  magnetic  field  quantities  will 
be  conserved  as  a  consequence  of  the  particle  equations  of  motion  and 
Maxwell's  equations.  However,  this  conservation  law  is  useless  in  a 

linearized  code,  which  does  not  calculate  nor  the  deviations  of  the 

electromagnetic  fields  at  second  order.  Thus,  we  have  taken  as  our  task 
the  derivation  of  conserved  quantities  which  are  sums  of  terms  bilinear  in 
the  first  order  quantities  computed  by  linearized  codes. 
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t. 


The  contents  of  this  article  are  organized  as  follows.  Section  2 
contains  a  derivation  of  the  linearized  perturbed  charge  and  current 
densities  for  general  coordinate  systems.  Specific  conservation  laws, 
e.g.,  continuity  of  charge  (or  mass)  and  a  Poynting  theorem  describing 
energy  flow,  are  derived  in  Section  3.  We  recover  the  results  of  earlier 
researchers6,7  in  the  cold  fluid  limit  (Sec.  4).  Illustrations  of  the 
conservation  laws  are  presented  in  Section  5. 


II.  PERTURBED  CHARGE  AND  CURRENT  DENSITIES 


In  this  section  we  derive  the  linearized  perturbed  charge  and 
current  densities.  We  shall  adopt  a  general  Klimontovich  representation 

g 

that  is  independent  of  the  particular  choice  of  coordinate  system.  For 
a  collisionless  plasma  satisfying  the  Vlasov  equation. 


if 

dt 


+ 


3v 


(1) 


where  f  =  f(x,v;t)  is  the  usual  phase-space  distribution  function,  the 
distribution  function  is  formally  given  by 


f(x,v,t) 


<5(x .  -  x)  5(7.  - 


where  the  sumnation  is  over  particles  and  x.(t),  v.(t)  satisfy 


(2) 


d_ 

dt 


d_ 

dt 


(3a) 


(3b) 
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The  brackets  In  Eq.  (2)  indicate  that  an  average  over  particle  initial 
conditions  has  been  performed,  which  is  equivalent  to  a  coarse-grained 
average  of  the  many-particle  distribution  function.  In  a  simulation 
plasma,  particles  are  finite  in  size;  and  the  shape  function  S( x ^  -  x) 

-►  -*■  Q 

replaces  the  delta  function  <5(x.  -  x).  Interpolation  of  charge  and 
current  densities  onto  a  spatial  grid  from  a  large  number  of  particles  per 
cell  and  solution  of  Maxwell's  equations  on  the  grid  for  the  macroscopic 

Q 

fields  accomplish  the  necessary  coarse-grained  average  in  the  simulation. 
Henceforth  in  our  analysis  we  drop  the  brackets  in  Eq.  (2)  and  presume 
that  the  averaging  is  understood. 

The  unperturbed  distribution  function  is  constructed  from  the 
unperturbed  trajectories, 

f(°)(x,v;t)  =  >T,  <s(x(0^  -  x)  -  v) 

i 

where  x1®)(t)  and  v^(t)  are  the  equilibrium  position  and 

velocity.  The  linearized  perturbed  distribution  function  is  then  given  by 


which  leads  to 


I  f  -  3°’)  5(;  -  'S0)) 

^  «(v  -  ^0))  -  ?S0)) 


(5b) 


upon  expansion  in  a  Taylor  series,  and  v(^  are  position  and 

velocity  displacements  relative  to  the  unperturbed  trajectories  (Fig.  1). 
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This  brings  in  a  Lagrangian  aspect  to  the  calculation  and  has  been  formally 
described  as  a  polarization  representation  in  the  excellent  treatment  given 
by  Bobroff.® 

Calculation  of  the  perturbed  charge  and  current  densities  from  Eq.  (5b) 
is  simple  and  direct.  The  linearized  perturbed  number  density  n^  for  each 
species  is 


„u)  ./d3;fm  ..£x*u> 

i 


(6) 


because  the  spatial  and  velocity-space  divergences  (3/3x)  •  and  (3/3v)  • 
of  the  particle  displacements  xl^,  x^  and  velocities  v^,  v^ 
are  identically  zero.  Similarly,  the  linearized  current  density  is 
given  by 


j(l)  a  q  /  yf 


« <»  -  -  e£  *1”  •  4 -  ;<°>j 


+  e 


(i) 


<5  l  X  -  X  : 


(0)' 


(7) 


Alternate  expressions  for  n^  and  can  be  obtained  with  use 
of  simple  vector  calculus  identities, 


.U) 


(8) 


We  observe  that  in  Eqs.  (8)  and  (9)  the  derivatives  (3/3x)  have  been 
coffinuted  with  the  suninatlon  over  particles.  This  allows  the  derivative  to 
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be  calculated  on  the  grid  in  the  simulation  code,  and  the  derivative  of  the 
shape  function  does  not  then  enter  explicitly. 

Use  of  the  identities  employed  in  deriving  Eqs.  (8)  and  (9)  yields  yet 
another  expression  for  the  linearized  current  density® 


;<”x  -<°)«( 


(10) 


Equation  (10)  has  the  important  computational  advantage  over  Eq.  (9)  in 
that  a  curl  of  a  vector  is  required  rather  than  a  divergence  of  a  tensor. 
Appropriate  care  must  be  exercised  in  evaluating  divergences,  curls,  and 
other  spatial  differential  operators  when  curvilinear  coordinates  are  used, 
which  in  our  personal  experience  has  been  a  frequent  source  of  difficulty 
in  both  analysis  and  code  development. 

To  evaluate  Eqs.  (6)  and  (7),  or  (8)  and  (9),  the  zero  and  first 
order  particle  trajectories  need  to  be  determined.  This  in  turn  requires 
the  electric  and  magnetic  fields,  which  typically  are  derived  by  solution 
of  Maxwell's  equations  on  an  Eulerian  mesh  given  sources  and  J^. 

The  forces  on  and  hence  the  acceleration  of  the  particles  are  deduced  from 
the  electric  and  magnetic  fields  evaluated  along  the  particle  trajectories. 

To  illustrate  the  polarization  representation  and  to  prepare  the 
groundwork  for  examples  that  are  considered  in  Sec.  5,  we  consider 
cylindrical  coordinates  for  a  system  that  is  axisyimtetric  at  zero  order 
(Fig.  1).  The  velocity  vector  is  given  by 


;  *  v(0)  +  v(1)  =  ;(°)  +  ^  7(1)  =  ^  (r0r0  +  z0z  +  Srr0  +  +  ^z) 


A,  •  A  .  *  • 


r0r0  +  ro0o9o  +  ZqZq  (^  “  ^g®g)**g  *  (Eg  +  Ey.0g)®g  +  ’  ( ^ ) 


tk 
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where  drQ/dt  3  909g,  deQ/dt  «  -90r0,  rQ  =  drQ/dt,  9g  =  d^/dt,  etc. 
From  Eq.  (TO)  identification  of  the  velocity  components  is  obvious. 


Ve 


+  W 


(12) 


(13) 


III.  CONSERVATION  LAWS 


A.  Continuity 

Integration  of  the  Vlasov  equation,  Eq.  (1),  over  velocity  space 
yields  the  continuity  equation 

#♦£•?■*  •  <"> 

Continuity  of  charge  (or  mass)  at  lowest  order  is  easily  demonstrated  with  use  of 

n^°)  3  <S(x  -  x(°>),  and  0^  *  e  v 6(x  -  xj°b  , 
i  i 

3n(0)/3t  3  -  (dxj°Vdt)  *  (3/9x)  6(x  ‘  xi^0)) 
i  ' 

3  -  ( 3/3x)  •  ^  ;;°)  /x  -  x|°A  3  -  (9/ax)  •  j(0)/e 
i  '  ' 

At  first  order,  differentiation  of  Eq.  (8)  with  respect  to  time. leads  to 


3t 


^°>\  -  3”  v*f°>  • 


.  05) 
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where  -  dx^/dt  and  vt°)  =  dlcj^/dt.  Because  (3/3x)  •  = 

-►  W  I  \ 

(3/ 3x)  •  xju  =0,  the  second  term  on  the  right  side  of  Eq.  (15)  can  be 
rewritten  as  -(3/3x)  •  vl^  xt^  <5(x  -  x^).  We  then  recognize 
that  Eq.  (15)  is  just  the  divergence  of  Eq.  (9)  multiplied  by  (-1),  i.e.. 


anil! 

at 


+ 


3_ 

-+■ 

ax 


(16) 


The  same  result  is  obtained  even  more  directly  by  using  Eq.  (10). 


8.  Poynting  Theorem 

Similar  conservation  laws  can  be  calculated  by  multiplying  the  Vlasov 
equation  by  increasing  powers  of  v  and  deducing  the  corresponding  moment 
equations.  Of  interest  here  is  the  energy  flow  in  a  linearized  simulation. 
The  linearized  Maxwell  equations  are 


(17a) 


7  X 


3^ 

c3t 


(17b) 


7  •  B{1)  =  0  (17c) 

From  the  scalar  products  of  Eq.  (17a)  with  E^  and  Eq.  (17b)  with 
we  obtain 

ft  (4^*4?)  ■  -  e(i)  •  j(”  .  as, 

This  is  a  Poynting  theorem  that  is  bilinear  in  perturbed  quantities. 

Because  Maxwell's  equations  are  "linear,"  the  same  manipulations  lead  to 
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Poynting  theorems  that  are  bilinear  In  field  and  current  perturbations  of 
arbitrary  order.  (In  typical  applications,  the  Poynting  theorem  constructed 
from  zero  order  quantities  may  be  especially  simple,  if  not  trivial,  due  to 
the  symmetry  and/or  simplicity  of  the  equilibrium,  e.g.,  8^  *  8QZ, 

Bq  =  constant,  and  E ^  3  3  0.) 

Monitoring  the  equality  of  the  left  and  right  sides  of  Eq.  (18)  in  a 
linearized  simulation  would  only  check  the  accuracy  of  the  numerical 
solution  of  Maxwell's  equations  and  would  give  no  information  on  the 
validity  and  self-consistency  of  the  particle  trajectories.  To  make  further 
progress  we  must  explicitly  evaluate  •  J^. 

The  linearized  equation  of  motion  is  given  by 


<£<i> 

dt 


:( 1 ) 
"tot 


:<°>  x  5(D 
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;<i>  x  l(0> 


(19a) 


where 
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The  scalar  product  of  Eq.  (19a)  with  vO)  gives 
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where  all  field  quantities  are  evaluated  at~x^.  Use  of  Eq.  (9)  for 


and  integration  by  parts  lead  to 


ft*  I<1>  •  j(’)  .  e  ■  l(1)  *  e  £:(°>  .  (*,  .  .  (2! 

Integration  of  Eq.  (18)  over  volume  and  substitution  of  Eqs.  (20)  and  (21)  yield 

Jf[/*  '*£■)*  E  i  ™»r)2]  *  kft ;  •  fc’)  *  *«m)  ■ 


Z[  *i,!  •  (*,”  •  $ )k°>  - .  (? . 

where  all  the  field  quantities  within  the  summation  over  particles  are 

-►f  Q\  ~ 

evaluated  at  xi  '  and  n  is  the  unit  normal  to  the  surface  containing 
the  volume  of  integration.  We  have  obtained  the  same  result  starting  from 
the  Low  Lagrangian. 10 

We  observe  that  Eq.  (22)  is  an  extended  Poynting  theorem  and  that  its 
right  side  is  not  zero  in  general.  In  fact,  the  right  side  cannot  vanish  if 
the  linearized  system  exhibits  growth  or  damping  and  appropriate 
boundary  conditions  are  adopted  so  that  the  Poynting  flux  (c/4n-)  x 
vanishes  at  the  surface  of  the  volume  (e.g.,  periodic  boundary  conditions,  a 
conducting  boundary,  or  vanishing  fields  at  infinity),  because  the  remaining 
terms  on  the  left  side  are  positive-definite.  The  right  side  of  Eq.  (22)  is 
bilinear  in  first  order  quantities  and  also  contains  zero  order  flow  and 
field  quantities.  This  allows  free  energy  at  zero  order  to  be  exchanged 
with  first  order  quantities. 
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C.  Energy  Conservation 

An  energy  conservation  law,  rather  than  a  Poynting  theorem,  can  in 
fact  be  derived.  Oberman  and  Dawson^  have  shown  that 


in  one  dimension.  Although  Eq.  (23)  is  often  described  as  an  energy 
conservation  law,  a  similar  equation  can  be  derived  from  invoking 
conservation  of  entropy,  as  can  be  demonstrated  by  considering  the  expansion 
of  the  quantity  -ftnf  to  second  order  in  small  perturbations  from  thermal 
equilibrium  and  by  employing  the  Vlasov  equation.  For  monotonic 
decreasing  f(v),  v"*3f^/3v  <  0  and  Eq.  (23)  measures  "energy"  transfer 
from  first  order  perturbation  fields  to  the  distribution  function.  For 
non-monotonic  distribution  functions,  i.e.,  v”1  3f^/3v  >  0  for  some  v, 
there  can  be  instability  and  the  energy  transfer  in  Eq.  (23)  is  in  the 
opposite  direction. 

It  is  important  to  note  that  f^  and  3f^/3v  are  not  usually 
calculated  in  a  linearized  particle  code,  and  their  evaluation  with  any  degree 
of  precision  is  difficult.  Numerical  calculation  of  (af^/Sv)-*  is 
especially  difficult.  In  the  case  that  3f  (°)4v  goes  through  zero  and 
changes  sign  as  a  function  of  v,  the  division  in  Eq.  (23)  may  not  be 
well-defined.  Of  course,  plasmas  with  non-monotonic  f^(v)  may  be  unstable 
and  are  therefore  of  particular  interest.  It  is  readily  seen  from  the  Vlasov 
equation,  however,  that  if  f^(x,v,0)  vanishes  at  v  3  v  where  3f(°V3v 
is  zero,  then  f^(x,vc,t)  =  0.  Thus,  it  is  plausible  that  Eq.  (23)  is 
well-behaved  for  non-monotonic  distribution  functions  provided  that  the 
initial  data  f^(x,v,0)  vanishes  at  the  zeros  of  af^/3v.  However, 
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whether  a  numerical  calculation  of  Eq.  (23)  in  a  simulation  would  be  well-behaved 
remains  a  serious  question. 

In  two-  or  three-dimensional  cases  when  the  equilibrium  distribution 
function  is  a  function  of  e  *  mv  /2  +  e$(x)  alone,  an  analogy  of  Eq.  (23) 
holds: 


H”)2 

23fi0,/3e 


If  f^/3e  vanishes  anywhere,  Eq.  (24)  is  then  subject  to  the  same 
limitations  on  the  initial  data  f^(x,v,0)  as  was  described  in  the 
preceding  paragraph. 

A  conserved  quantity  can  also  be  derived  from  a  Lagrangian  formalism 

based  on  an  extension  of  the  Low  Lagrangian.10  In  a  paper  by  Galloway  and 
13 

Kim  ,  the  Low  Lagrangian  for  a  warm  plasma  and  finite-amplitude 
electromagnetic  waves  is  derived  and  generalized.  The  single-particle 
Lagrangian  is  given  by 


</  -  \  mv2  -  q(4>  -  |  •  A)  ,  (25) 

where  v  *  dx/dt  and  $  and  A  are  the  electrostatic  and  vector  potentials, 
respectively.  The  field  Lagrangian  density  is 

/  *  h  if*  *  it)'  4  x  *) 2]  •  (26) 

At  this  point  Galloway  and  Kim  introduced  an  Eulerian  phase  space 
description  in  preparation  for  averaging  of  the  Lagrangians.  To  distinguish 
Eulerian  phase  space  variables  from  the  polarization  variables  used 
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elsewhere  in  this  paper  we  define  Galloway  and  Kim's  variables  by  (X,V)  * 

(Xg  +  X^,  Vq  +  Xj).  Here  Xg  and  Vg  are  reference  position  and 

velocity.  X^  =  X-|(Xg,Vg,t)  represents  a  small  displacement  of  a 

phase  space  element  that  in  the  absence  of  perturbations  would  be  located  at 

Xg  and  have  velocity  Vg  (see  Fig.  2  in  Ref.  13).  X^  =  [(3/3t)  +  Vg  •  (3/3Xg)  + 

-►  -*■  + 
aQ  •  (3/3Vg)]X1  is  the  velocity  perturbation,  with  aQ  the  particle  acceleration 

in  the  unperturbed  fields.  Xg  and  Vg  are  independent  fluid  variables,  and  X^ 

and  X^  are  dependent  fluid  variables.  When  (Xg,  Vg)  =  (x^,  v^) 

the  functional  values  of  X^  and  X^  coincide  with  the  values  of  the 

-W  1  )  -+{ 1  \ 

polarization  variables  xv  '  and  vlA/  introduced  earlier. 

Galloway  and  Kim  obtained  an  approximate  Lagrangian  by  expanding  the 

Lagrangians  and  .2^  in  perturbation  expansions  with  respect  to 
-►  -►  •+ 

perturbations,  X^,  X^,  <j>,  and  A^.  The  particle  Lagrangian  is  summed  over 
plasma  particles  by  assuming  that  the  plasma  is  distributed  in  an  approximately 

P 

continuous  fashion  and  by  converting  the  sum  ofi?  over  particles  to  an 
integral  f d3Xd3Vf(X, V,t]2>P  where  f(X,V,t)  is  the  phase  space  number 
density.  Because  the  number  of  plasma  particles  in  a  cell  is  conserved  under 
perturbation. 


d3X  d3V  f(X,V,t)  «  d3Xgd3V0fg(Xg,Vg) 
where  fg(Xg.Vg)  is  the  equilibrium  phase  space  number  density.  Neither 

^  -f 

fg(Xg,Vg)-  nor  the  equilibrium  fields  Eg  and  Bg  is  time  dependent.  The 
system  Lagrangian  at  second  order  in  perturbation  quantities  is  equal  to 

L2«  f ^3xgd3Vgfg(Xg,Vg)^2  +  f  2  •  (27a) 


where 
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2  z  *  \  m  Xx2-  q|(X1  *  V)<(,1  +  \  (X^:  W)*Q  -  JJ  •  *  V)^  +  £  (X^  :  W) 

-  r  •  [*i  *  <*i  •  Sj  27b) 


(27c) 


Introduction  of  the  action,  SJ  dt  L^t),  and  setting  its 

variation  equal  to  zero,  6^  =  0,  lead  to  Euler-Lagrange  equations  that 

are  the  equations  of  motion  for  the  first  order  plasma  displacement  and 

13 

electromagnetic  fields: 


The  variations  <SX^,  5<^  and  vanish  at  the  endpoints  t  =  tj  and  t^, 

but  are  otherwise  independent  and  arbitrary  functions.  Hence,  to  satisfy 

6 Ig  a  0  each  of  the  square  brackets  in  Eq.  (28)  must  vanish  separately,  which 

13 

then  yields  the  Euler-Lagrange  equations. 


1 
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We  have  derived  a  conservation  law  from  the  Galloway  and  Kim 
Lagrangian  by  using  the  fact  that  L2  in  Eq.  (27)  has  no  explicit  time 
dependence.  We  consider  variations  induced  by  an  infinitesimal  time 
translation,  t  -»•  t  +  At,  so  that  6X1  =  At(3X1/3t)  [the  variations  need 
not  vanish  at  t  =  t^  and  t2  in  Eq.  (28)].  If  the  unperturbed  orbits  and 
fields  are  constrained  to  satisfy  the  Euler-Lagrange  equations,  then 


rtz  rl  2  dL0  r  1 

/  dt6L2=At/  dt^ 

=  At 

|L2  ( t2 

)  - 

L2(tl)J 

ll 

r  0  0  •+ 

3X: 

3^P 

O 

t2 

At  I  dJX0dJVQf0  (XQ,V0)  - 

"Tt  * 

tl 

9 

which  upon  rearranging  terms  allows  us  to  identify  an  energy-like  conserved 


quantity  E2, 


'<Wo<VV" 


X:  -  L2  =  constant. 


(29) 


where  L2  is  given  by  Eq.  (27).  It  is  important  to  note  that  E2  is  not 
the  second  order  system  energy  and  that  E2  is  not  composed  of  only 
positive-definite  quantities. 

It  is  illuminating  to  consider  an  example  of  the  conservation  law,  Eq. 
(29).  For  Langmuir  waves  in  a  one-dimensional  homogeneous  plasma  with 
^(x)  =  <j>jexp(-iwt  +  ikX)  +  c.c.. 


from  which  it  follows  that 
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me(co  -  kVQ)‘ 


Z-L  exp(.-io3t  +  ikXg)  +  c.c. 


X1  +  V0  ^)X1 


me(w  -  kVQ)‘ 


exp(-  icot  +  ikXg)  +  c.c. 


The  disDersion  relation  follows  from  Poisson^s  equation , 


±*£  fdV  W  -W  C 

me  J  0  (w  -  kVn)2  mek  J 


o  JJTZ  kV0) 


For  a)  =  cjr  +  1y,  the  real  and  imaginary  parts  of  the  dispersion  relation  yield 


i  -  4’e‘  f  it  r  It  )  <“r  *  kVo)2^2 
1  •  IT-  /  dVo(vo>  77 v - 

6  J  ("r  -  k*o)  *  »* 


(32a) 


r  2y(u  -  kVQ) 

o-J^oWu  '°~ 

(“r  -  kVo)  * 


(32b) 


We  substitute  the  results  of  Eq.  (30)  in  Eq.  (29)  and  average  over  a  spatial 
wavelength  X  =  Z^/k  to  obtain 


/c  Y  I*  \2a2n  fdV  f  (V  )  2UV(aV  ~  kV 

Y2/  H7  ^i1  J  dWV  ( ZTT^fT? 


This  is  proportional  to  the  imaginary  part  of  the  dispersion  relation  Eq.  (32b) 
and  therefore  vanishes.  Thus,  despite  the  possible  growth  or  damping  of  the 
Langmuir  wave,  is  zero  and  is  indeed  conserved. 
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Equation  (29)  differs  markedly  in  philosophy  from  our  Poynting 

theorem,  Eq.  (22).  In  the  Galloway  and  Kim  approach  the  plasma  is  conceived 

to  be  a  continuous  Vlasov  fluid,  while  our  computational  plasma  is  an 

ensemble  of  discrete  particles.  Calculation  of  fg(Xg,Vg)  and  3X^  /8t 

in  a  linearized  particle  simulation  would  require  a  phase-space  grid.  While 

the  sum  of  the  particles  at  all  velocities  within  a  spatial  grid  cell  is 

large  enough  to  insure  adequate  statistics  in  particle  codes,  if  phase  space 

were  additionally  subdivided  into  velocity  bins  there  would  typically  be 

very  few  particles  in  each  cell.  The  fg(Xg,VQ)  calculated  in  the 

simulation  would  then  only  be  approximately  time  independent  because  of  the 

significant  statistical  fluctuations  in  the  number  of  particles  wandering  in 

>  ^ 

and  out  of  the  phase-space  cell.  The  quantity  ^/at  or  [VQ  •  (3/3XQ)  + 

■+  ■*  -*■  -*■ 

aQ  •  (3/3Vq)]X2  would  require  finite  differences  of  X1  determined  locally 

in  phase  space  and  would  also  be  difficult  to  calculate  accurately  because 

of  statistical  problems.  It  is  therefore  our  conclusion  that  Eq.  (29)  is 

an  interesting  and  valuable  analytical  result,  but  it  does  not  lead  (at 

present)  to  a  practical  particle  code  diagnostic.  Equation  (23),  its 

generalization  Eq.  (24),  and  Eq.  (29)  are  much  more  practical  as 

diagnostics  in  simulations  in  which  the  linearized  Vlasov  or  drift-kinetic 

equation  is  solved  directly  for  f^(x,v,t). 

Although  Eq.  (22)  does  not  provide  a  conserved  quantity,  as  do  Eqs. 

(23),  (24),  and  (29),  it  possesses  certain  essential  advantages  that  lead  us 

to  recottmend  its  use  as  a  diagnostic  in  a  particle  code.  Equation  (22),  or 

its  time  integral,  contains  quantities  that  are  readily  available  in  a 

linearized  particle  code  and  is  easily  calculable  with  good  precision. 

Furthermore,  the  breakdown  of  the  right  side  of  Eq.  (22)  into  three  terms 
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affords  valuable  diagnostic  information  on  the  energy  flow  in  and  out  of  the 
first-order  quantities  and  checks  the  accuracy  of  the  solution  of  both 
Maxwell's  equations  and  the  particle  equations  of  motion. 

IV.  COLD  FLUID  LIMIT 

In  the  cold  fluid  limit,  all  particles  within  the  immediate 

*+■  •* 

neighborhood  of  a  field  position  x  are  presumed  to  share  the  same  x 
vi  ' ,  and  vi  ' .  At  a  finite  distance  away,  these  values  may  differ 
due  to  finite  spatial  gradients,  but  they  are  independent  of  particle 
index.  Then  it  follows  that  the  linearized  perturbed  number  density 
becomes 

„d) . .  (3/3i)  s(: .  ;(»>) . .  (3/3;j .  . ;(«>> 

i  ' 

-  -  (3/3x)  •  (nQ^(1))  ,  (33a) 


and  the  current  density  is  similarly  given  by 
J(1)/e  *  nQv(1)  -  (3/3x)  •  (n^1^0)) 

=  ( 3/3 t)  (nQx(1))  +  (3/3x)  x  (nQx(1)xv(0))  .  (33b) 

These  expressions  agree  with  those  found  in  Refs.  6  and  7.  Bobroff, 

it 

Haus,  and  Kluver  pursued  this  type  of  analysis  further  and  derived  a 
small-signal  power  theorem  of  the  form  (3/3t)W  +  (3/3x)  •  S  s  0  for  a  cold 
continuum  fluid,  Eqs.  (4.11)  and  II. 7)  of  Ref.  7.  Our  Eq.  (22)  and  the 
energy  conservation  law,  Eq.  (29),  obtained  from  the  Galloway  and  Kim 
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formalism  can  be  shown  to  be  mathematically  equivalent  to  the  power  theorem 

"  7 

of  8obroff,  Haus,  and  Kluver  in  the  cold  fluid  limit. 

14 

In  subsequent  papers  Sobroff  examined  energy  conservation  for 
plasma  models  based  on  the  Vlasov  equation  and  on  amodel  of  a  system 
composed  of  many  cold  beams.  The  Vlasov  model  gave  the  result  presented  in 
Eq.  (23),  but  the  beam  model  is  rather  different  and  interesting.  Consider 
N  beams  with  number  densities  and  velocities  n^,  u.,  i  *  1,  2,  ...N  and 
an  electric  field  E  each  of  which  satisfy 


3ni  + 

ST*7’  <Vi>  =0 


^Ui 

sr  *  (“i  •  7K  ■  <S>. 


N 

f  +  47  £  eiVi  * 0 

1=1 

If  this  system  is  linearized  about  an  equilibrium  n^),  uj®),  E^)  and  the 
equations  for  the  first  order  perturbations  nP),  uP),  E^1)  are 
deduced,  then  there  results  the  energy  conservation  law 


3_ 

3t 


*ui» 


(34) 


2  "iP 

i  X 


(1)7  X  up)  -  nP)  7  xu^ 

i  i  i  i 


Thus,  for  a  one-dimensional  system  the  right  side  vanishes  and  a  true 

14 

conservation  law  is  obtained,  as  was  calculated  by  Bobroff.  It  is  also 


evident  from  this  generalization  of  Bobroff's  analysis  to  beams  in  two-  or 
three-dimensional  systems  that  there  is  a  conserved  quantity  only  for 
irrotational  flows;  otherwise,  Eq.  (34)  is  in  the  form  of  a  Poynting  theorem. 


V.  EXAMPLE:  MAGNETIZED  ELECTROSTATIC  HYBRID  CODE  IN  CYLINDRICAL  POLAR  COORDINATES 


In  this  section  we  consider  an  example  that  illustrates  the 
calculations  of  the  preceding  sections.  Byers15  has  developed  a 
linearized  electrostatic  code  which  solves  the  Poisson  equation. 


,  i  ,V’> 


.  .  £  tU)  .  i  .  £  *»>  .  ,„e  („jl>  .  ,  ,35) 

and  advances  the  trajectories  of  particle  ions  self-consistently  in  an 
applied  magnetic  field  B  *  BgZ.  The  first-order  electrostatic  potential, 
the  derived  electric  fields,  and  the  charge  desnities  are  assumed  to  vary 
as  <P  8  <T(r,t)exp  ime  +  c.c.  where  $(r,t)  is  a  complex  amplitude. 

The  electron  fluid  response  is  given  in  the  low-frequency 
approximation  (3/3t  <  <  eBg/mgc)  with  the  electron  polarization  drift 
assumed  negligible. 


3n<°> 


(36) 


Because  Eq.  (36)  gives  an  expression  for  3rTe/3t,  Byers'  code  in  fact 
determines  the  electrostatic  potential  by  using  the  time-derivative  of 
Poisson's  equation,  Eq.  (35).  Derivatives  with  respect  to  r  are  represented 
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on  a  uniform  radial  grid  as 


— 7 
I3r‘ 


bd) 


Jl)  ,,(1)  .  ,(1) 

2Ar2  ~ 


g(1)  , 

r,J 


1_*0) 

3r  y 


Jj 


*0)  .  *(D 
!i±i.  r 

2ar 


and 


Ed) 

fc9,j 


3  ^(1)  3  .  jm  (1) 
J  J 


n36 


where  j  is  the  spatial  grid  index. 

The  equilibrium  is  axisymmetric  and  trivial  in  the  B  *  0  limit  (0  « 
plasma  pressure/magnetic  pressure):  E^  *  0,  BQ  is  uniform  and 
constant,  and  the  ions  execute  simple  cyclotron  motion  in  the  (r,9)  plane. 
In  the  presence  of  electric  field  perturbations,  the  zero  and  first  order 
trajectories  are  described  by  the  solutions  of  Eqs.  (10)  and  (11).  In  this 
case  z  is  a  completely  ignorable  coordinate,  and  we  set  z^  =  ?z  =  = 

=  0  for  all  ions.  The  ion  equations  of  motion  are 


d_ 

dt 


a  V(0) 

m^c 


x  B 


(37) 


at  lowest  order  and 


(38) 
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at  first  order.  The  right  sides  of  Eqs.  (37)  and  (38)  are  evaluated  at  the 
lowest  order  ion  position  x^.  This  requires  interpolation  from  the 
radial  grid  to  the  particle,  for  which  linear  cloud-in-cell  interpolation  is 

a 

used.  The  particle  equations  of  motion  are  integrated  using  a  centered 
3  g 

leap-frog  scheme.  * 

An  essential  simplification  results  from  the  axi symmetry  of  the 
equilibrium  and  the  exp(im0)  variation  of  the  first  order  electric  field. 

All  ions  at  the  same  radius  and  same  velocity  vector  (vr,vQ)  have  the  same 
zero  order  trajectory  except  for  time- independent  displacements  AS.  These 
in  turn  only  lead  to  multiplicative  factors  exp(imA6)  in  the  first  order 
velocities  and  positions.  Therefore,  a  complete  angular  distribution  of 
particles  need  not  be  tracked.  It  is  sufficient  to  follow  one  particle  per 
Ztt  radians  if  comp lex -valued  amplitudes  are  used  to  represent  all  first 
order  quantities.  However,  it  is  convenient  here  in  our  formal 
manipulations  to  replace  the  single  super-particle  by  a  continuous  angular 
distribution  of  point  particles  at  each  radius. 

The  ion  number  density  perturbation  is 


,0)  = 


*  -  (a/ax) 


2X 


ro  +  5eVi 


S(r 


r(0)j 

r-1 - 6(0-9|O)),  (39) 


where  (3/3x)  •  *  [r  r“*( 3/3r)r  +  0r"^(3/30)]*.  The  grid  representation 
S(r  -  rj°)')  replaces  <S(r  -  r(0^),  3/30  =  im,  and  (3/3r)  is  evaluated  on 
the  grid  in  the  code.  The  summation  on  particles  in  Eq.  (39)  is  performed 
to  generate  a  grid  quantity  before  the  divergence  is  calculated.  The 
first  order  current  density  is  determined  in  Eq.  (13), 
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i 


where  (3/3x)x  =  [r(3/3r)  +  9r“*(3/39)]x.  The  expressions  for  n^  and 
given  in  Eqs.  (39)  and  (40)  obviously  satisfy  the  continuity  equation. 
The  Poynting  theorem,  Eq.  (22),  is  simple  for  this  system. 


where  v^  3  (vj^rg  +  4°^0o^i  and 


(1) 


(1) 

9 


+  z 


9  r39 


E0)  iaE0) 

r  re 


+  ^9  r30 


+  ^9 


(41) 


(42a) 


(42b) 


■*■(01  -*■  -*■ 
Because  electron  inertia  is  neglected,  E'  '  3  0,  and  only  the  electron  E  x  B 

response  is  retained,  there  is  no  electron  contribution  to  Eq.  (41),  i.e., 

E  *  Jg  3  0.  Furthermore,  for  this  system  it  is  easily  shown  from  Eqs.  (7) 

and  (20)  that  the  ion  contribution  to  the  right  side  of  Eq.  (41)  is 

equivalent  to  -  f  d^x  E^1^  •  [  -  e  ^vj1)  <5(8  -  <5(r  -  r^)/rj, 

l 

which  is  more  easily  computed  in  the  code. 
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Figures  2,  3,  4,  and  5  illustrate  the  consistency  of  linearized 
simulations  with  the  continuity  equation  (16)  and  the  Poynting  theorem  (41) 
for  both  stable  and  unstable  plasmas.  The  first  example  of  a  stable  plasma 
(Fig.  2)  consists  of  a  kinetic  equilibrium  with  spatially  overlapped  ion 
orbits.  The  subsequent  examples  (Figs.  3,  4,  and  5)  use  a  plasma 
corresponding  more  to  a  fluid  equilibrium  with  all  particle  orbits 
encircling  r  *  0  and  satisfying  conditions  for  a  rigid  rotation.  In  a 
stable  plasma  the  amplitudes  of  the  various  normal  modes  are  constant  in 
time  and  are  determined  by  the  initial  conditions.  Because  there  is  apt  to 
be  a  multiplicity  of  frequencies  for  a  given  wavenumber,  especially  for  an 
energetic  plasma  in  a  magnetic  field,  we  should  not  be  surprised  to  observe 
the  beating  of  stable  normal  modes  that  is  evident  in  our  energy  diagnostic 
(Fig.  2).  In  the  unstable  case  (Figs.  4  and  5),  the  kinetic  and  electric 
field  perturbations  grow  exponentially  in  time  driven  by  the  free  energy 
made  available  by  the  density  gradient  and  the  ion  distribution  function, 
which  is  monoenergetic  at  any  radius.  The  Poynting  theorem  diagnostic 
provides  clear  evidence  of  the  coupling  of  the  velolcities  v(0^  to  the 
differential  change  in  the  first  order  electric  field  associated  with  the 
first  order  ion  displacement,  x^  •  (3/3x)E^  at  x  *  x(°). 

We  believe  that  the  continuity  equation  and  the  Poynting  theorem 
provide  essential  diagnostic  checks  on  the  accuracy  of  the  linearized 
simulations  considered  here  and  also  furnish  valuable  insight  into 
microscopic  processes.  Furthermore,  we  have  found  that  these  conservation 
laws  have  played  a  vital  role  in  debugging  and  developing  new  linearized 
codes  and  in  isolating,  subtle  difficulties  associated  with  nonuniform 
equilibria  and  curvilinear  coordinates  in  particular. 


“cit/2ir 

Figure  2.  Linearized  simulation  of  a  stable  plasma  with  azimuthal  mode  number  m  *  2, 
2  2 

pi'/^ci  *  **  an<^  dens1ty  scale  length  Rp  *  4a.  The  perpendicular 

velocity  distribution  is  monoenergetic  vi  *  aiU)ci’  and  the  ion  orblts 

overlap  in  space.  There  is  a  conducting  wall  at  r/a^  *  2.  (a)  E  = 

f  d?x  E^  /8ir,  J  =  f  dt'/*d3x  E^  •  and  the  total  energy 

o  2  *  2 

T  =  E  +  J  vs  time,  (b)  E  =/ d3x  E(1)  /Sir.  <  =  J  (l/2)m1vj1)  , 

0-  =^dt-S^0)  G,  ■  0/31)  E<'l]  ■jrSt'/’  d3x  !">  •  Cj(') 

-  e  £  5(9  -  6(r  -  r(®h/r],  and  the  total  energy  T  = 

1  / 

E  +  K  +  0  vs  time. 


■ 


} 


Figure  3.  Linearized  simulation  of  a  stable  configuration  of  axis-encircling 

2  2 

Ions  with  azimuthal  mode  number  m  *  2,  3  1  at  r  *  0, 

and  density  profile  n^  =  oq  exp(-r^/rg).  The  zero  order 
trajectories  satisfy  v^  *  -r^  where  3  eB/m^c. 

There  is  a  conducting  wall  at  r  3  2rg.  (a)  E,  J,  and  the  total  energy  T  = 
E  +  J  vs  time,  (b)  E,  K,  j',  and  the  total  energy  T  =  E  +  K  ♦  J'  vs  time. 


Figure  4.  Linearized  simulation  of  an  unstable  (Imuj*  0.2  wci)  configuration  of 

|t  2  2 

axis -encircling  ions  with  *  1.5,  and  otherwise  the 

same  parameters  and  circumstances  as  in  Fig.  1.  (a)  E,  J,  and  T  *  E  + 


Figure  5.  For  the  simulation  described  in  Fig.  4,  the  real  and  imaginary  parts  of 
the  complex  amplitudes  for  SnPVst  and  -(3/3x)  •  vsr  r/rQ. 

at  a  particular  time.  The  equivalence  of  these  quantities  to  high  degree 
demonstrates  conservation  of  charge  and  mass. 
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F.  DIGITAL  FILTERING  IN  TIME  AND  HIGH  FREQUENCY  DAMPING 

Dr.  W.  M.  Fawley  (LLL)  (Prof.  C.  K.  Birdsall) 

Crystal  et  al.  (1979)  have  pointed  out  that  numerical  dissipation 
may  be  inserted  into  plasma  simulation  codes  by  backward  biasing  of  the 
field  equations  and/or  the  equations  of  motions.  If  two  quantities  are 
related  by  the  expression,  (e.g.,  velocity  and  electric  field) 


si  *  c,<0 


Crystal  et  al.  suggest  that  replacement  of  the  usual,  centered  equation 


„n+T  _  _n- j 


M  oi  -At 


,n+i  _  _n-+ 


/  3  n+1  1  n-K 

( Tp  +  Vs  ) 


will  lead  to  frequency-dependent  numerical  damping.  Their  algorithm  has 
the  especially  nice  property  of  the  imaginary  part  of  the  computed  complex 
frequency  scaling  quadra t i ca 1 ly  with  the  normal  mode  frequencies  (e.g., 
electron  plasma  frequency)  of  the  system.  Thus,  one  might  introduce 
this  dissipation  in  order  to  damp  the  high  frequencies  in  the  system  (e.g., 
ope^  while  keeping  low  frequency  behavior  intact,  and  then  take  much  larger 
time  steps  than  normal  stability  requirements  would  allow  (i.e.,  ui  At'2). 

However,  the  above  algorithm  requires  implicit  knowledge  of  quanti¬ 
ties  at  times  n+l  in  order  to  advance  other  quantities  to  time  n++.  Crystal 
et  al.  suggest  implicitly  solving  for  these  advanced  quantities  via  matrix 
inversion.  While  in  principle  one  might  attempt  an  iterative  solution  for 
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an+'  in  Eq.  (3),  our  experience  has  been  that,  as  the  time  step  increases 

toward  w  At  *  2  where  m  is  the  maximum  normal  mode  frequency  of  the 
max  max 

system,  the  number  of  iterations  necessary  for  a  stable  solution  approaches 
infinity.  Thus,  we  agree  with  Crystal  et  al.  that  an  implicit  solution 
must  be  used  for  Eq.  (3)  —  a  point  made  by  A.  B.  Langdon  numerous  times. 

Unfortunately,  if  an+'  depends  upon  individual  particle  properties 
such  as  position,  the  necessary  matrix  inversion  for  implicit  solution  be¬ 
comes  impracticable.  One  variation  of  the  backward-biasing  present  in  Eq. 
(3)  appealed  to  us  which  would  minimize  the  need  for  iteration.  Rather 
than  use  as  large  a  time  separation  as  2At  on  the  right  side  of  Eq.  (3), 
why  not  decrease  this  separation  to  25At  with  5  «  1  as  in  the  following  set 
of  equations: 


where  0<e<1  is  the  biasing  factor. 

The  hope  in  using  Eq.  (4)  is  that  since  vn  ^  and  En  are  known,  then 
one  can  make  a  very  accurate  prediction  for  En±l5  for  6  small.  We  employed 
the  following  algorithm: 


.n-i 


■ 


(7) 


V 


n 


n-± 


iAt 


a  e"-4 

m 


(8) 


n±6 

x 


n 


x 


±  <5«At*vn 


(9) 
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The  algorithm  from  time  (n-$)  to  (n+i)  is  shown  in  Fig.  1.  The  sequence  is: 

„(n-l)*  n-i  n  _n  _n-£  n  n±5  _n*  .  .  , 

E  -*■  v  -*■  x  -»■  E  -*•  E  -"v  x  -+•  E  .  The  equations  are  time-centered 

except  for  expression  (9)  which  deviates  to  first  order  in  <$. 

If  one  analyses  this  algorithm  by  assuming  all  quantities  vary  as 

exp(-icut)  and  that  the  normal  mode  frequency  is  to  ,  one  obtains  for  to  At  <  1 

o  o 

with  £*0.5,  5*0.2,  A t  =  1.0  (see  Figs.  2S3), 


■j  5 

to  "  to  ^  0. 04to  At 

r  o  o 

2 

to.  £  .06(o  At 
o 

i.e.,  the  method  produces  frequency-dependent  damping.  Figures  4  and  5 
show  the  effect  on  <or  and  to.  of  varying  5  while  keeping  all  other  para¬ 
meters  constant. 

Introducing  the  algorithm  into  ESI  was  not  difficult  and  the 
only  significant  additional  memory  requirements  were  supplying 

the  array  space  for  En”  .  For  to  At  <2.0,  the  algorithm  worked  as  expected 

o 

in  that  high  frequencies  were  damped  much  more  than  low  frequencies.  We 
also  experimented  with  arbitrarily  turning  the  damping  on  and  off.  One  ex¬ 
ample  is  shown  in  Fig.  6  where  the  damping  was  "on"  for  0<t<  150  and  "off" 
for  t>150.  As  expected,  the  electrostatic  energy  in  mode  32  exponentially 
damped  for  t<  150  and  then  remained  constant  for  greater  times. 


i-lh 

xn*jL 


Vn_2 


vn  xr 


'  n*i 
Vn  2 


FIG.  1  Leapfrog  integrator  with  micro  (5At)  and  normal  (At)  time  stepping, 

showing  detail  about  time  tn.  The  micro  steps  are  used  to  alter- 

n* 

nate  the  high  frequency  content  up  the  electric  field,  with  E  used 
to  advance  the  velocity  from  time  n-i  to  n+i. 


i 
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FIG.  2.  The  real  part  «  of  the  dispersion  relation  is  plotted  versus 
the  driving  frequency  (e.g.,  u>  ) .  Here  At  =  1.0  sec,  e-0.5, 
and  5*0.2.  The  straight  1  i  ne  w  *  u)  is  also  shown. 
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IG.  3  The  imaginary  plot  w.  is  plotted  versus  u >  for  the  same  para¬ 
meters  as  in  Fig,  2.  Note  the  quadratic  increase  injcu.jwith 
increasing  u>  which  leads  to  preferential  damping  of  high 
frequencies . 
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Unfortunately,  the  algorithm  is  not  stable  for  u>  At  > 2.0.  The 

o 

problem  is  probably  in  expression  (7)  where  En  should  be  replaced  by  its 

biased  counterpart,  which  again  requires  an  implicit  solution.  It  is 

possible  that  iteration  may  be  more  useful  in  this  algorithm  than  that  of 

Crystal  et  al.  To  work  around  this  instability  problem,  we  experimented 

with  the  following  "hybrid"  time-step  scheme:  For  each  long  time  step 

( co  At,  >2)  taken  by  the  code  which  excites  high  frequency  oscillations, 

o  i  ong 

take  many  (Nshort)  time  steps  using  the  algor i thm  above  to  damp  out 

these  unwanted  oscillations.  Normal  modes  with  oj«  2At , '  should  remain 

long 

unaffected  by  the  varied  time  step.  In  practice,  we  discovered  the  fol¬ 
lowing  behavior.  Each  long  time  step  leads  to  an  exponent i a  1  growth  in 
the  electrostatic  energy  contained  in  high  frequency  modes  and  the  argu¬ 
ment  of  the  exponential  scales  1 i nearly  wi th  At^^  (see  Fig.  7).  Thus, 
the  ratio  of  Nshort^shor  t  to  N^At^  must  be  equal  to  or  greater  than 

some  constant  for  stable  behavior  in  all  modes.  For  j  *1.0,  At  .  *  1 , 

o  snort 

£*0.5,  5*0.2,  N  .At  .  /N,  At,  S  20. 

short  short  long  long 

That  is,  in  using  this  hybrid  scheme  one  gains  less  than  5%  in 

the  total  time  simulated  for  a  given  number  of  time  steps  as  compared  with 

a  normal  simulation  just  using  At  ,  .  Therefore,  unless  a  solution  is 

found  for  the  a)QAt>2.0  instability,  the  damping  algorithm  does  not  appear 

to  be  useful  for  very  long  time  step  simulations  despite  its  simplicity  of 

coding  and  frequency-dependent  damping  properties. 

It  is  to  be  noted  that  we  still  have  gained  some  advantage  in 

being  able  to  work  up  to  a  ]tvl,5,  with  damping  increasing  quadratic- 

ally  with  a)  At,  producing  the  desired  (undamped)  physics  at  «u)  .  The 

pe  pe 

ultimate  goal  of  being  able  to  allow  m  At  > 2  still  eludes  us. 

pe 
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FIG.  7-  The  jump  in  the  l°9ar '  thm  0,:  che  electrostatic  energy  in 

a  given  mode  from  taking  a  single  time  step  with  ^t]ong  ~^0 
is  plotted  versus  At,  •  The  approximate  linear  relationship 
indicates  that  the  number  of  short-time  steps  employing  time- 
filter  i  ng  requ  i  red  to  damp  out  this  jump  varies  directlv  with 
the  length  of  the  jump  in  it. 


Crystal,  T.  L.  ,  Oenavit,  J.,  and  Rathman,  C.  £. ,  Comments  on  Plasma 
Phvsics  and  Controlled  Fusion.  5,  17  (1979). 
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Section  II 

CODE  DEVELOPMENT  and  MAINTENANCE 


A.  ESI  CODE 

No  special  report  this  quarter. 

B.  EMI  CODE 


No  special  report  this  quarter. 

C.  EZOHAR  CODE 

No  special  report  this  quarter. 

D.  RINGHYBRID  CODE 

Alex  Friedman 

Several  errors  were  made  in  the  original  implementations  of  the 
code;  these  have  not  affected  results  for  ion  rings  significantly,  since 
they  are  small  for  currents  at  radii  large  compared  to  the  cell  size,  in 
general.  For  plasmas  with  current  near  the  axis,  such  as  FRM  plasmas,  it 
may  be  important  to  use  the  correct  expressions  (derived  elsewhere  in  this 
report) . 

The  first  error  was  a  failure  to  note  the  vectorial  nature  of  the 
term  (g^ •  ^7) ;  due  to  the  azimuthal  derivatives  of  the  unit  vectors,  two 
terms  result: 


a  contribution  to  J 


1 

rk 


£9k9k 


a  contribution  to  j],  * -e v  ? / r? 

9k  9k  rk  k 
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The  next  error  was  a  failure  to  recognize  the  radial  variation  of 
v°  for  each  superparticle,  necessary  if  there  is  to  be  no  shear  to  the  flow. 
The  first  consequence  of  this  is  the  appearance  of  an  extra  term  in  the 
first-order  current: 


a  contribution  to  J  .  «  -s  .  9.  . 

0 k  rk  k 


The  second  consequence  is  that  azimuthal  current  densities  J°  j'  must  be 

0  0 

computed  in  a  different  manner  than  n°: 


J,:  divide  by  2 nr  (particle)  rather  than  2irr  (cell). 

0  O 


J  :  divide  by  2irr  (particle)  rather  than  2irr  (cell) 

0  o 


These  changes  were  implemented  into  the  code,  and  V-J  now  closely 
.  1 

approximates  -n.  Results  regarding  stability  of  ion  layers  and  rings  were 
not  altered  greatly. 


E. 


RADIAL  SIMULATION  CODE  ES1RB 


No  special  report  this  quarter 

F.  RADIAL  CODE  NOTES  (R . Re , R2 , ReZ) 

1.  R0  Area  Weighting  of  Charge  (Method  A) 
C .  K.  S i rdsa  1  I 


Charges  are  to  be  weighted  to  an  R0  grid.  Let  them  be  rod  shaped 
clouds  with  curvilinear  shape  Ar  by  rAP  uniform  in  z  as  shown  in  Fig.  1. 
The  grid  has  cells  Ar  by  r A<|> .  The  cell  in  which  the  particle  appears  is 
given  inner  radius  r.,  outer  radius  r.  ,  and  angular  grids  i>  ,i.  .. 

i  i+l  a  j  j+1 

The  coordinates  of  the  particle  r  , p  are  taken  to  be  the  arith- 

P  P 

metic  mean.  The  curvilinear  cloud  is  taken  to  have  charge  q  coulombs  /  m, 

P 

with  uniform  charge  density  at  the  instant  of  weighting.  As  the  area  of 


at  the  next  instant;  no  account  is  made  for 

(expands  jn  angu|ar  lenqth. 

(contracts 

In  an  axisymmetric  system,  with  no  p  variation  there  are  just 
charged  shells  of  thickness  Ar,  as  in  Fig.  2.  Ignoring  the  grid,  the  force 
on  shell  2  of  density  q„,  due  to  shell  1  producing  field  E^  is  (see  Fig.  3). 


the  cloud  i 

is  (r2 

-r2  )  A 

outer  inner 

(decreases 

as  r 

I  i  ncreases 
(decreases 

\i  ncreases 

P 

the  energy 

(loss 

as  the  cloud 

l,ga  i  n 

F  (on  2  due  to  1 ) 
r 


h 


E.  dV 


(E  is  rad i a  1  ,  , 

p  is  in  Cou lombs/nv) 


where 


o 


(q  is  in  Coulombs/m) 


L 


! 

i 

! 

i 
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0 


FIG.  1. 


Cloud  of  charge  q,  a  rod  parallel  to  z,  uniform 
located  at  r.#p(*),  dr,  6$  from  grid  point  ij. 
are  weighted  to  nearest  grid  points.  The  cloud 
Ar  is  curvilinear,  with  A$,  Ar  fixed.  The  grid 
A<t>,  Ar,  fixed. 


in  density, 
Shaded  areas 
shape  rA<J>  by 
uses  the  same 


Cloud  shell  of  charge  q,  a  cylindrical  shell  centered  on  origin  0, 

uniform  in  density,  located  at  (*)r  >  Sr  from  grid  point  r. .  The 

shell  is  the  rod  of  Fig.  1,  with  angular  width  A<j>  made  to  be  2*; 

the  radial  width  is  t  =  ir.  For  r  *Ar/2,  the  inner  shell  radius 

P 

is  at  ra0;  for  r  <Ar/2,  the  shell  inner  radius  moves  outward, 

P 

crossing  the  outer  radius  for  rp*0;  that  is,  the  shell  folds 


through  the  origin. 
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So 


F 

r 


2tte  r 
o 


2irr  drdz 


q1P2(b  -a) 


E 

O 


dz 


for  p2  independent  of  r.  (As  the  particle  moves,  with,  say,  shell  thick¬ 
ness  kept  constant,  then  within  a  shell  p.  changes;  the  assumption  here  is 
that,  at  any  given  position,  p  is  independent  of  r.)  This  expression  is 
awkward  to  use  as  p^  changes  as  the  particle  moves;  this  form  can  be  im¬ 
proved  by  changing  p^  into  a  charge/length  quantity  which  does  not  change 
as  the  particle  moves,  writing 


/,  2  2. 

P2ir(b  -  a  ) 


P22irt 


m 


p,2irtr 

2  P 


Hence, 


F^  (on  q2  due  to  q1  inside  of  q2) 


qiq2 

2irVP 


(Newtons/m) 


=  q2Er(at  r  due  to  q^  ; 


that  is,  use  Er  as  measured  at  the  center  of  the  charge  being  moved.  Of 

course,  E^  may  be  due  to  many  q.  all  lying  within  shell  q2;  simply  replace 

q .  above  with  See  Fig.  4  for  thin  shells 

i“1  1 

There  is  also  a  self-force  in  that  a  cylindrical  shell  tends  to 

expand,  in  contrast  to  a  planar  sheet,  which  does  not  (with  some  exceptions, 

like  the  sheet  in  question  not  midway  between  sheets  of  opposite  sign,  or 

at  t*,  etc.)  That  is,  a  shell  at  r  produces  a  field  between  itself  and 

P 

charge  of  opposite  sign  at  r>rp’  suc^  t*iat  there  is  a  force  on  the  shell 

of  value  /  pEdV;  this  force  for  p  constant  over  the  shell  (E  is 
•4he  1 1  r 


Variation  of  with  neutral  sets  of  thin  charged  shells.  The 
jump  in  D  is  (Cou iombs/m^)  which  is  q/2irr .  Note  the  con¬ 
trast  with  planar  sheet  models,  with  the  same  jump  at  each  sheet 
and  no  variation  between  sheets. 
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(Newtons/m) 


The  factor  in  parentheses  varies  from  1/3  to  1/2  as  a/t  varies  from  0  to  « 
(from  rod  on  axis,  to  thin  shell).  Hence,  the  self-force  varies  from  1/3 
to  1/2  that  of  one  shell  of  charge  q  inside  the  shell  of  interest  of  charge 
q;  for  many  charges  inside  a  shell,  the  self-force  may  be  negligible  (un¬ 
less  the  net  charge  enclosed  is  comparable  to  q) . 

2.  Method  B 

Weighting  in  R6  may  be  done  using  the  recipe  implied  in  Fig.  1. 

This  will  be  done  straightforwardly.  Then  the  limit  A0-*2tt  will  be  made, 

producing  ring  (or  shell)  charges  and  radial  weighting  to  be  compared  with 

Method  A  (see  QPR  III  for  this  year). 

The  total  charge  to  be  weighted  to  the  nearest  four  grid  points 

is  q  (as  q  is  a  rod,  the  units  are  Coulombs/meter).  The  area  of  an  element 

of  q  of  angular  width  a  extending  between  r.  and  r  is 

33  3  inner  outer 


;(area)  =  f fr  drd$  =»  ^-(r2  -  r?) a  =*  (r  )(width)a 

JJ  2  o  i  mean 


where  r  *  (r.  +  r  )/2  and  (width)  ■  r  -  r. . 
mean  i  o  o  i 


The  fractional  area  to  be  assigned  to  a  grid  point  is 
S(area) 


(r  ) (width)a 
mean 


total  area 


r  ArA$ 
P 
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Using  the  symbols  in  Fig.  1,  6r»r  -r.,  <J>  - <j> . ,  we  may  construct  r 

pi  P  J  mean 

and  width  for  the  four  grid  points  and  hence,  the  weights.  This  is  done 
in  Table  1. 

Adding  all  of  the  [(r  )(width)a]'s  does  indeed  turn  out  to 

mea  n 

be  rpArA<)>,  a  first  check  (total  area  or  charge  conserved). 

For  rings  (or  shells),  let  Ad-*-2ir.  Then  the  weight  (a  fraction 
of  charge)  assigned  to  r.  is  q., 

qi  *  q  (ri  +T')  (Ar_5r)/rpAr 

and  that  assigned  to  r .  ,  is  q.i1f 
*  i+l  i+l 

qi+1  *  q  (ri+  Al~2~  r  )  (<Sr)/rpAr  . 

We  next  compare  these  weights  with  Method  A  of  QPR  111,  this  year. 

Table  1 
R-e  A kzol  Wzighti 


poi  nt 

r  inner 

rou ter 

r 

mean 

width 

a 

i,j 

BB 

BbH 

EH 

Ad>-6<t> 

i  +  1  ,  j 

HEM 

V*r+¥ 

EE3 

B 

A  $-6$ 

i+l , j+1 

B 

B 

ar 

5$ 

i  ,  j  +  1 

B 

mm 

B 

Ar-6r 

5d> 

■it 

5r»r  -r .  and  64  *  <)>  -4 . 
pi  P  j 


-  76  - 


Using  the  subscripting  of  Fig.  1,  with  grids  at  r.,r.+^  and  particle 
located  at  r  between  r.,r.^,, 

p  i  i  +1 


p-6r  “i 


the  results  of  Method  A  are 

q .  »  q (r . ) (Ar-6r)/rpAr 

q!+1  “  d(r.+))(^r)/rpAr  . 

The  differences  are  slight.  Method  B  has  the  first  factors  as  the  mean  of 
(r.,r  )  or  (r  ,r..,)  where  Method  A  has  r.,r.  Method  B  is  derived  for- 

i  p  p  i  +  1  It-*-! 

mally,  is  taken  to  be  more  accurate,  and  is  recommended  as  a  charge  in  ex¬ 
isting  codes.  Be  careful;  note  whether  charge  or  densitv  is  being  weighted. 

*****7**  *******  ********  ********  ******** 

These  notes  are  being  included  somewhat  randomly,  to  be  collected 
later  in  a  more  ordered  fashion,  with  tests,  criticisms,  recommendations,  etc 

*******  *******  ********  ********  ******** 

Previous  QPR,  Fig.  1  on  Page  9^  is  perhaps  more  clearly  understood 
if  the  ordinate  q(r.)  is  replaced  by  qr(r.);  that  is,  r.  is  being  varied, 

J  j  ' 

not  r..  The  particle  at  r.  is  moving  past  the  fixed  grid  point  r.. 


RJET  DEVELOPMENT 
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No  special 


report  this  quarter 


H.  SOFTWARE  DEVELOPMENT  (RUN .  SOLVER) 


RUN 


H.  Stephen  Au-Yeunq  and  Niels  Otani 


Eemq  the  first  -an try  of  a  LIB  file.  RUN  executes  any 
program  within  the  LIS  file.  The  purpose  of  having  RUN  is  to 
encourage  jsers  to  out  their  programs  into  a  library.  RUN  is 
:n  the  CRAY- 1  computer  and  can  be  obtained  by  typing: 
rfilem  read  1222  .Cray  runtESCI  end  t  v 


Th i 5  document  corresponds  to  the  RUN  version  of  FcDrary  11, 
1?S0.  _dter  versions  of  RUN  will  be  stored  in  “ ILSN  directory  .cray 
of  user  numoer  1222.  The  user  should  per  iod icai 1 y  check  the  date  of 
these  -'lies  ifilem  now  1222  .cray  run'  to  see  if  the  program  has  been 
updated,  ""he  file  RUN  is  a  LIB  file:  it  contains  the  latest  source 
as  well  as  the  latest  version  of  this  document  a t : on .  To  obtain  this 
document,  type; 

lib  runtESCx  run.  doc  tESC'end  t  v 

i- e tout  Cusc]  run  doc  C  turn.  3  box  bnn  run  t  v 


Suppose  RUN  i s  the  f  i-st  entry  of 
cor<a  ms  proarams  '.ike  CCMOUT  ana  FREY, 
may  execu'e  70MOUT  and  TREX  as  follows: 
•dll  10  ! 


he  LIE  'lie  ALL.  which 
As  an  example,  the  user 


run  user 
comout  user 
comou t 
run  user 
f  rex 


-ccmout  ucb  endwith.  s  s.  box  o22  source 
Nend 

Ycomoutk  :  file-id  is  comout  -  rxuebka 
-frex  f ICSroOx  6  t  SO  box  022  fr^Ofile 
Processing  file:  f 10Sro0x 


run  user 


end 


Note  that 
‘he  premot 
. o  pub  1  i  c 


ALl.  will 
' -  ' .  If 
tile,  for 


not  be  terminated  unt 
the  oregram  entered 
example!.  RUN  will  a 


1  "end”  is  entered  after 
s  not  in  the  1  ibrary 
so  execute  it. 
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SCL-.ER 

■-i.  Svepnen  'u-'isunq 


SOLvER  has  been  updated  on  Narcn  is.  '.5S0.  “ha  new  version 
hots  (ne  following  cncnctes: 

v  1 '  5CL-ER  ncu  works  on  tne  nsu  LTE  for ~"od*.s.  "or  mesa  mo  would 
:kj  to  ■.•so  mo  old  lDR  may  dot  mo  o.d  version  of  SOLVER  by 
1  yp  mg  -  -filom  rood  1020  .classic  so  1 . or ,S5u ' end  ‘  v 

or  mo  CRPY. 

-O'  S  row  variable.  'MRk .  has  boon  addod  to  mo  NAMELIST  read 

so  that  ‘.ho  user  may  specify  tne  marker  for  tr.e  curve  vs)  of 
nu  nor  roots  or  function  plot*,  slid. 


NV'Rk  «  ■*■'.  «s  points  connected  and  d  symbol  everu  ith  point. 

-i  *N  points  not  connected  and  a  sumeol  at  every  ith 
OOint . 

0  *s  points  connected  with  no  symbols  drawn 
voef oul 1 1  . 

Sromer  row  variable.  FNREhL.  nas  been  ddded  to  the  NAMELIST 
recta  .o  or.aole  SOLVER  to  1  ock  for  only  pure  roal  roots. 

FNRE.SL  *  1  »'  to  force  oil  appro-,  imat  ions  to  oil  the  roots  to 
bo  rod  1 . 

0  *'  allow  complev  -'oo  t  s  .default' . 
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bility  of  Field  Reversed  Ion  Rings". 
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I  .  PURPOSE 

This  conference  was  held  in  order  to  facilitate  an  exchange  of 
ideas  and  results  on  plasma  particle  and  particle-fluid  hybrid  codes  in 
an  informal,  unpublished  manner,  with  frank  and  open  discussion  among 
about  fifty  people.  The  only  printed  record  is  the  collection  of  copies 
of  transparencies  used  by  the  thirty-five  speakers.  In  keeping  with  this 
informality,  the  transparency  copies  were  distributed  to  the  participants 
only  and  are  not  to  be  quoted  without  explicit  permission  from  the  author. 

In  December  197**  we  organized  a  similar  unpublished  conference  in 
Berkeley  at  the  request  of  Robert  Price  of  ERDA.  The  present  conference  is 
a  follow-up  and  is  due  in  large  part  to  the  interest  of  David  Nelson  (Chief, 
Fusion  Theory  and  Computer  Services  Branch,  Div.  Applied  Plasma  Physics, 
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Office  of  Fusion  Energy,  00E).  Dr.  Nelson  presented  his  view  from  DOE  of 
contributions  by  computation  and  simulation  toward  furthering  plasma  re¬ 
search  leading  to  fusion  reactors;  he  received  considerable  feedback  from 
participants.  We  are  most  grateful  to  him  for  his  active  and  stimulating 
part i ci pat  ion. 

The  conference  is  also  due  to  the  pressure  of  accumulated  interest 
among  those  who  work  with  particle,  and  hybrid  codes,  who  wished  to  exchange 
new  ideas,  methods  and  results,  in  an  informal  setting.  By  no  means  was 
this  meeting  intended  to  conflict  with  the  Ninth  Conference  on  Numerical 
Simulation  of  Plasmas  scheduled  for  June  30  -July  2,  1 980  at  Northwestern 
University,  sponsored  by  Professors  J.  Denavit  and  G.  Knorr,  which  covers 
a  much  wider  area  and  is  to  be  published. 

Lastly,  the  particle-hybrid  code  community  had  concern  that  the 
very  real  contributions  of  particle  and  particle-fluid  simualtions  to  the 
understanding  of  fusion  plasmas,  in  support  of  both  theory  and  experiment, 
might  just  have  been  underestimated  in  Washington.  This  conference  is,  in 
part,  a  reaction  to  recommendations  of  the  DOE  committee  for  computer  time 
allocations  for  FY80  which  included:  major  labs  receiving  63%  of  their  re¬ 
quests,  but  universities  receiving  only  40%;  the  decision  that  large  part¬ 
icle  pushing  or  kinetic  codes  can  be  afforded  only  sparingly;  basic  plasma 
theory  was  given  approximately  a  20%  absolute  cut;  ability  to  follow  out 
unexpected  results  was  cut  to  a  bare  minimum;  studies  of  alternate  con¬ 
cepts  were  cut  more  than  were  mainline  studies.  The  challenges  to  the  com¬ 
munity  are  clear:  make  our  contributions  known;  make  our  codes  more  effi¬ 
cient  in  terms  of  useful  physics  per  unit  of  computer  time  (optimize,  use 
minimum  number  of  dimensions  and  particles,  etc.);  check  code  physics  by 
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analysis  (time  and  spatial  grid  effects,  fluctuations)  and  by  comparison 
with  known  linear  and  nonlinear  results  so  as  to  increase  confidence  in 
the  simulations. 

Editorially  we  note  that  the  field  has  matured  considerably  from 
the  early  desk  calculator  work  of  Buneman  and  Hartree  on  magnetrons  in 
the  early  1 9^*0 * s ,  the  one  dimensional  electrostatic  plasma  work  on  modern 
computers  by  Buneman  and  Dawson  in  the  late  1950‘s,  to  the  current  2d  and 
3d  fully  electromagnetic  codes  of  Langdon,  Buneman  and  others.  The  early 
science  has  become  more  of  an  art.  Simulation  for  magnetic  fusion  has  now 
become  strongly  applications  oriented,  with  DOE  devoting  the  bulk  of  FY80 
computing  resources  to  the  direct  support  of  existing  experiments  and  the 
design  of  next  generation  devices,  including  calculations  of  transport, 
impurities,  heating,  stability,  equilibrium,  coil  design,  etc.,  a  mile¬ 
stone  indeed.  Coupled  with  this  new  emphasis  is  the  current  plateau  in 
large  computer  time  for  magnetic  fusion,  with  no  relief  planned  until  at 
least  summer  1 98 1 .  These  factors  present  a  challenge  to  all  of  us  to  be 
more  efficient  and  more  effective  in  our  simulations. 


Co-on.cUncLtosiA : 


C.  K.  (Wed)  Set d&ati  A£ex  F fLczdman 

U.  C.  BeAkeZzy,  17ec.embe/t  17,  1 979 
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II  . 


Final  Schedule  for 

Informal  Conference  on  Particle  and  Hybrid  Codes  for  Fusion 
December  10-11,  1979,  Napa,  California 


Conference  coordinated  by  Charles  K.  (Ned)  Birdsall,  Alex  Friedman, 
assisted  by  Singer  Fletcher,  at  the  Electrical  Engineering  and 
Computer  Science  Department,  University  of  California,  Berkeley  94720 

There  will  be  four  main  sessions.  Talks  will  be  20,  15  and  10 

minutes  with  adequate  time  for  discussion.  If  certain  subjects  draw 
great  interest,  additional  time  will  be  available. 

Of  general  interest  is  the  efficiency  of  particle  and  particle- 
fluid  codes,  in  terms  of  physics  output  per  unit  of  computer  resource 
(time,  memory,  volume  of  output,  etc.).  Of  comparable  interest  is  the 
progress  toward  working  at  lower  and  lower  frequencies,  with  larger  mass 
ratios,  etc. 


MONDAY,  DECEMBER  10 

Session  IA  -  8:30  am  to  12:00  noon  Chairperson:  C.  Nielson 

*  Long-time-averaging  (LTA).  Particle  simulation  of  slow  transport 
phenomena . 

*  Particle  MHD  vis  a  vis  fluid  MHD  and  other  fluid  codes. 

*  Large  time  step  problems  (mi/me>>1,  omega  dt>>1,  digital  filtering, 
stiff  equation  integrators) . 

gafiSE*  Pag 

1  J.  Denavit,  T.  L.  Crystal,  C.  E.  Rathmann,  J.  L.  Vomvoridis,  "Applica-  ~i 

“  tions  of  Long-Time-Scale  Particle  Simulations." 

2  B.  Cohen,  R.  P.  Freis,  T.  A.  Brengle,  "An  Orbit-Averaged  Particle  Code."  22 

3_  W.  Fawley,  "Low  Pass  Filtering  in  Time."  44 

4  D.  Anderson,  "Toward  a  Full  3d  Hybrid  Transport  Code."  51 


5  T.  Tajima,  J.  M.  LeBoef,  F.  Brunei,  J.  M.  Dawson,  "Recent  Efforts  in  56 


Particle  MHD  Code  Development." 

7  H.  Okuda,  "Steady  State  Drift  Turbulence  and  Anomalous  Diffusion,  3d."  101 

6_  C.  K.  Chu ,  "Vortex  Ring  Formation  by  Vortex-in-Cell  Method."  93 

9_  C.  Tull,  "Code  Exchange  Mechanics."  121 

^  3.  McNamara,  ''Remarks  on  the  Use  of  Integrals  of  Motion"  j_2Q 


Paper  numbers  here  correspond  to  paper  numbers  in  the  proceedings'  book,  which 
was  sent  to  all  participants. 


Session  IB  -  1:30  pm  to  5:00  pm  Chairperson:  C.  K.  Bird  sail 


*  David  Nelson;  view  from  DOE. 

*  Improvements  in  simulation  1974  to  1979.  and  future. 

*  Code  production  vs  development  running  times;  allocations  and  priori¬ 
ties. 

*  Hardware,  e.g.  array  processors,  graphics,  class  VII  computers. 

*  Efficient  3d  simulation;  3d  grid  vs  2d+Fourier  representation. 


C.  K.  Birdsall,  Comments  on  1974  Berkeley  Meeting;  hopes  for  this  meet¬ 

ing. 

D.  Nelson,  "Role  of  Particle  and  Hybrid  Codes  Present  and  Future;  Com¬ 

puter  Availability;  Possibility  of  Adding  Special  Purpose  Computers 
for  Particle  Codes." 

A.  B.  Langdon,  "Tradeoffs  Among  Code  Development  vs  Hardware  Costs  vs 
_  _ _ Elapsed  Calendar  Time;  Future  Hardware  Needs;  ZED  Postprocessor." 


Page 
~ I4T 

145 


J.  Kulp,  "High  Performance  Array  Processor  for  LIST  Machine;  .  __ 

Architecture,  Impact  on  Particle  Simulation." 

B.  Moore,  W.  Drummond,  "Particle  Simulation  o  n  the  VAP."  159 

T.  3rengle,  N.  Maron ,  3.  Sutherland,  "Use  of  an  Array  Processor  165 

With  PD P-10. ” 

R.  3erman,  "Macrocell  Algorithm  for  Efficient  Particle  Pushing, 

for  CRAY  and  AP  in  2d,  2§d,  3d."  iba 


■(T.“TT"(5ieng,  H.  Okuda,  "3d  Simulation  of  Trapped  Electron  Instabilities  181 

in  Toroidal  Systems." 

0.  Buneman ,  "Data  Management  for  a  Million-mode  3-d  e-m  Code;  Use  of  196 

Tetrahedral  Mesh." 

T.  Tumolillo ,  "MEEC-3D:  Description  of  an  Existing  Self  Consistent  Par-  210 

tide  Pusher." 
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TUESDAY,  DECEMBER  11 

Session  IIA  -  8:30  am  to  12:00  noon  Chairperson:  J.  Denavit 

*  Linearized  codes. 

•Modified  particle  codes  (use  of  linear  susceptibility  or  Boltzmann 
response.  Id  stretched  to  2d  or  3d). 

*  Quasineutral,  hybrid,  and  Darwin  codes  (applications  to  confinement, 
compression,  equilibrium,  stability,  and  transport). 

*  Buildup  and  plasma  trapping,  neutral  beam  injection,  and  wave  heating 
(laser-pellet  plasmas,  pinches,  mirrors,  tokamaks). 

Page 

W.  W.  Lee,  H.  Okuda,  "Particle  Simulation  Models  for  Low  Frequency  233 
Microinstabilities, 

A.  G.  Sgro,  "Hybrid  Simulation  of  Non-MHD  Phenomena."  251 

D.  Hewett,  "A  Global  Method  of  Solving  the  Electron  Field  Equations  in  a  263 
Zero-Inertia  Electron  Hybrid  Simulation." 

22  D.  Winske,  "Particle  Simulation  of  Reversed  Field  Configurations."  281 

23  R.  Mason,  "Monte  Carlo  (Hybrid)  Model  for  Electron  Transport  in  Laser  288 

Plasmas." 


24  J.  Byers,  "Id  Linearized  Particle  Model  for  Tandem  Mirrors  and  Field  301 

™  Reversed  Mirrors." 

25  A.  Friedman,  J.  Denavit,  R.  N  Sudan,  "Linearized  3d  Hybrid  Simulations;  320 

Ergodic  Orbits  in  Simulation." 

26  V.Decyk,  "Diagnostics  for  Bounded  Plasma,  with  Applications."  350 

27  B.  Cohen,  N.  Mar on ,  G.  R.  Smith,  W.  M.  Nevins,  "DCLC  Simulations  with  a  365 

Stretched  Id  Code." 

28  R.  Huff,  "Particle  Hybrid  Codes  on  the  CHI  Computer."  339 


Session  IIB  -  1:30  pm  to  4:00  pm  Chairperson:  0.  Buneman 

•  Inhomogeneous  plasmas  (fluctuations,  initialization  techniques  in  2d 
and  3d) . 

•  Undesirable  instabilities  in  warm  plasma  simulations  (e.g.  multi-beam 
and  multi-ring)  . 

•  Grid  effects  (e.g.  curvilinear  coordinators,  div  B  nonzero). 


20 

30 


W.  Nevins,  "Fluctuations  in  Inhomogeneous  Systems." 


V.  Thomas,  C.  K.  Birdsall,  "Alias  Growth  in  Hybrid  Oscillations  due  to  4^9 
Initiation  at  k  A  x  *  ir." 


31 


A.  Drobot,  A.  Palevsky,  "E-M  Simulation  of  Strongly  Radiating  Systems. 
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4:00  pm  to  5:00  pm  Moderator:  B.  Cohen 

*  Panel  Discussion:  When  to  use  a  particle,  fluid,  or  hybrid  code  (on 
none  at  all).  All  participants  invited  to  contribute. 

Panelist:  J.  Denavit,  A.  B.  Langdon,  B.  McNamara,  C.  Nielson,  H.  Okuda 


B.  Godfrey,  "Electromagnetic  Numerical  Instabilities  in  Two-Dimensional 
Relativistic  Beam  Simulations." 


A.  Sternlieb,  "Coupling  of  Particle  Codes  to  Electric  Circuits." 

J.  Poukey,  J.  P.  Quintenz,  "Preliminary  Simulations  of  Ion  Beam  Neutral¬ 
ization." 

Y.  Chen,  "Multi-Beam  Instability  Interference  with  Lower  Hybrid  Drift 
Instability." 
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III.  PANEL  DISCUSSION  * 

(8.  Cohen,  moderator:  J.  Denavit,  A.  B.  Langdon,  B.  McNamara, 

C.  Nielson,  H.  Okuda) 

Plasma  simulation  is  one  tool  in  the  attack  on  understanding  plas¬ 
mas  which  will  eventually  result  in  the  design  of  fusion  reactors.  Simula¬ 
tion  using  thousands  of  particles  vies  for  usefulness  with  simulation  using 
a  few  fluids.  Particle  codes  can  deliver  the  full  kinetic  behavior,  which 
is  sometimes  needed;  fluid  codes  use  parameters  from  theory,  particle 
simulation,  and/or  experiment,  and  deliver  average  or  long-term  information. 
Particle-fluid  hybrids  attempt  to  use  the  best  features  of  both,  mixing 
fast  and  slow  time  scales. 

Particle  codes,  in  some  quarters,  have  gained  a  reputation  for 
using  more  computer  time  per  unit  of  useful  physics  than  used  by 
fluid  codes.  This  reputation  is  considered  undeserved 

by  those  particle  simulators  who  work  with  limited  budgets,  use  the  least 
number  of  dimensions  for  the  problem  at  hand,  use  the  least  number  of  par¬ 
ticles  necessary,  and  use  hybrid  simulation  where  applicable.  However, 
there  appears  to  be  room  for  further  optimization.  Those  with  larger  com¬ 
puting  budgets,  please  follow. 

Okuda  addressed  the  progress  made  in  understanding  anomalous  trans¬ 
port  indigenous  to  tokamaks,  moving  from  an  early  fully  dynamic  (hence,  ex¬ 
pensive)  model  to  guiding  center  electron  models  in  3d,  both  electrostatic 
and  magnetostatic,  with  long  time  steps  (less  expensive),  toward  a  3d  tor¬ 
oidal  model  for  steady  state,  including  turbulence  and  transport,  poloidal 
divertors,  and  hybrid  heating.  These  models  are  particularly  useful  for 

This  section  was  written  from  my  notes  and  has  been  checked  over  by  the 
panelists  (C.  K.  Birdsall). 
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low  3  plasmas.  He  also  gave  an  example  in  which  the  low  frequency  ion 
acoustic  instability  can  drive  the  high  frequency  two-stream  instability 
unstable  by  nonlinear  effects;  this  is  a  warning  that  the  long-time-step, 
time-averaged  low  frequency  method  may  break  down  for  certain  plasma  para¬ 
meters.  Whether  or  not  this  kind  of  process  can  occur  for  other  low  fre¬ 
quency  instabilities  is  not  known.  In  any  case,  one  must  be  extremely 
cautious  when  modeling  low  frequency  phenomena  by  discarding  high  fre¬ 
quency  responses  since  they  can  make  a  large  change  in  what  one  is  looking 
for. 

Denavit  observed  that  hybrid  codes,  3d  particle  codes,  with  time 
filtering,  long  time  scales,  all  appear  to  be  problem  dependent.  In  con¬ 
flict  with  this,  there  is  the  need  for  optimization  via  assembly  language, 
which  implies  loss  of  flexibility.  We  are  fortunate,  in  particle  simula¬ 
tion,  to  be  able  to  work  on  interesting  problems  and  to  obtain  useful  re¬ 
sults.  Simulation  identifies  physical  mechanisms,  suggests  theoretical 
work  and  rejects  ideas  that  do  not  fit. 

Nielsen  stated  that  simulation  must  continue  in  order  to  keep 
theorists  honest;  simulation  is  the  only  economical  means  of  checking 
theory  directly  open  to  them.  The  question  is  not  whether,  but  how  to  con¬ 
tinue,  effectively.  We  must  pay  more  attention  to  macroscopic  work;  as 
we  do,  we  are  likely  to  find  that  distinctions  between  particle  and  fluid 
work  will  disappear.  The  new  hardware  is  revolutionary,  but  there  still 
is  an  open  question  as  to  using  one  big  computer  or  many  small  machines; 
his  view  is  that  special  computers  are  still  too  small,  too  complicated  to 
use. 

McNamara  presented  a  list  of  unresolved  problems  as  follows: 


Tokamatis 


1.  Anomalous  electron  transport 
—  is  it  really  classical? 

2.  Nonlinear  effects  of  ballooning  and  tearing 

Field  i?Cl ’&Ual 

1.  Non-axi symmetr i c  states  in  RFP 

2.  Equilibria  and  stability  of  large  orbit  FRM 
Tandem  MLixoxa 

1.  Nonlinear  saturation  of  cyclotron  modes 

2.  MHO  stability  and  ballooning 

3.  Energy  transport  in  tandems,  thermal  barriers,  etc. 

His  upbeat  view  is  as  follows: 

I.  Simulation  is  very  effective  and  increasingly  realistic. 

Therefore,  push  major  design  and  physics  efforts  for  fusion 
systems . 

Provide  more  software  and  graphics  support  for  designated 
activities.  Others  benefit  by  fallout. 

3.  Encourage  more  collaboration  —  from  three  day  advisory 

groups  for  new  problems  to  joint  efforts  and  publications. 

■f.  Improve  the  numerical  analysis  input  —  without  a  S200K  tax 
to  pay  for  i  t. 

3.  Write  more  reviews  on  the  state  of  the  o-t t  and  advertise 
good  results  more  vigorously. 

3.  Long  term  goals  —  Design  a  fusion  reactor:  Are  we  getting 
the  re? 

McNamara  also  presented  a  chart,  computational  attack  on  plasma  parameters, 
as  attached.  Note  the  marriages  of  all  kinds,  across  fuzzy  boundaries,  par¬ 
ticle  plus  fluid;  micro  plus  macro;  dynamics  plus  equilibria,  etc. 

McNamara  ended  by  repeating  the  need  for  more  advertising  of  accomplishments. 
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THE  COMPUTATIONAL  ATTACK  ON  PLASMA  PARAMETERS 
(B.  McNamara) 
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Langdon  noted  that  the  MFE  computational  effect  has  grown  greatly 
in  the  last  five  vears.  He  repeated  his  advice  of  the  1 9 7 ^  meeting  that, 
if  we  want  to  work  with  long  time  scales,  implicit  time  integrations  will 
be  required  which  mav  not  be  attractive,  but  are  needed ,  because  explicit 
methods  will  end  up  unstable  at  large  jAt.  He  suggests  turning  to  con¬ 
tinuum  methods  for  modeling  kinetic  physics  -  painful,  but  needed  (Rod 
Mason  has  done  some,  with  multigroup  transport).  He  is  interested  in  par¬ 
tially  kinetic  models. 

Cohen  concluded  the  oane  I  discussion  with  the  observations  that 
scaled  variables  will  probablv  alwavs  be  with  us,  that  having  computer 
variaoles  match  experiment  is  oroDaolv  impossible  in  the  foreseeable  future. 
He  presented  a  summary  cnart  on  hvprid  s i mu  1  a t i ons .  attached.  His  work 
and  that  of  others,  has  moved  particle  simulation  into  phvsics  on  trans¬ 
port  time  scales. 

SirdsaM  sketched  on  the  blackboard  the  DOE  hierarchy  of  models, 
with  increasing  expense  treading  left  to  right). 
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making  the  point  that  these  models  interact,  are  mutually  supportive,  and 
are  seldom  wholly  independent;  over-emphas i : i ng  one  tool  or  removing  an¬ 
other  might  well  reduce  the  efectiveness  of  the  remaining  svstem. 
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Kinetic  detai 1  requires  particles  [or  f(v)  description] 

a  / R  ■  @(1±)  large-orbit  equilibrium  &  transport, 

sp  e.g.,  field  reversal;  a  =v  /u> 

3  s  is  cs 

a  /X  ■  @(1±)  microinstability,  e.g.,  drift  waves 

S  W3VC 

Xn/L  a  (f)(1±)  sheaths 

u  n 

u  -  kj  V|  + 1(  a)  .>  +  puibounce  *  ®  Landau  resonance,  e.g.,  wave  (±)  dis- 

^  01  ou  .  u  sipation,  trapping,  overlap  of  resonance 

When  kinetic  detail  is  unnecessary ,  elimination  of  particles  and  fine  time 
&  space  scales  and  adoption  of  fluid  description  can  greatly  extend  compu¬ 
tational  power.  me/n'j  *1  often  separates  ions  &  electrons. 

u)/o)  ,  a  /R  «  1  drift  approximation 

cs  s  p 

ou/tucs  ,  ag/Rp  >:>  *  unmagnetized,  straight-line  orbits 

<v>  «  io/k  hydrostatic 

<v>  »  tu/k  adiabatic 

linear  phenomena  linear  dielectric,  xs  (u>  ^  i  3 1 ,  k) 

Payoff:  Analytic  reduction  of  fluid  description,  appropriate  to  specific 
parameter  regimes  and  problems,  (l)  result  in  large  savings  of 
computer  memory  &  time,  and  (2)  greatly  extend  the  range  of  para¬ 
meters  accessible  to  simulation. 


Pi rections : 

1.  Simulation  of  transport  with  hybrid  orbit-averaged  code 

2.  Hybrid  simulation  of  microinstabilities  with  more  realistic 
parameters  and  geometries. 
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IV.  CONCLUSIONS 

(from  C.  K.  Birdsall's  notes) 

Many  of  these  are  already  in  the  Pane)  Discussion. 

The  197^  conference,  for  the  same  purposes,  covered  a  number  of 
topics  which  are  now  pretty  well  resolved,  e.g.,  particle  weighting  battles 
are  mostly  past,  with  linear  and  quadratic  weighting  widely  used  and  under¬ 
stood  (multipole  use  isolated);  hybrid  codes  were  just  starting  then  and 
now  exist  in  a  variety  of  forms;  3d  codes,  then  done  by  one  or  two  groups 
only  are  now  done  by  several;  one  special  purpose  computer  was  mentioned 
then,  with  many  now. 

However,  several  authors  from  1974  gave  up-dates  in  1979  on  pro¬ 
blems  still  not  fully  resolved,  such  as  linearized  particle  codes,  long¬ 
time-scale  particle  simulations,  the  plea  for  use  of  implicit  schemes  for 
large  time  steps,  and  3d  economical  codes.  The  1974  conference  called  for 
more  use  of  hybrids  (happening),  more  emphasis  on  simulations  including  in¬ 
homogeneous  density,  magnetic  field,  temperature,  etc.  (happening)  and 
for  honest  use  of  two  species,  large  m./m^  (happening  somewhat).  Perhaps 
the  long-term  needed  to  resolve  these  problem  areas  indicates  their  tough¬ 
ness  and  rasies  the  question  of  pursuing  them  further  -  some  of  them  obviously 
must  be. 

David  Nelson's  view  from  DOE,  paper  #10,  was  provocative.  The  loud 
and  clear  messages  were  that  magnetic  fusion  computing  capability  is  fixed, 
with  no  relief  expected  before  1981,  with  increasing  allocations  to  device 
specifics  and  less  to  large  particle  pushing  or  kinetic  codes.  His  list  of 
problem  areas  needing  more  attention  from  particle  simulators  was: 
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Oivertors,  axi- + non-axi-symmetric  sheaths,  potentials,  boundary 
layers,  ambi polarity, 

Plasma  edge  including  shadow  of  limiter,  neutrals,  atomic  physics, 
reflux,  metal  surface, 

Transport  in  non-axi symmetric  systems, 

RF  heating  including  nonlinear  effects 

Optimization  of  gyrotron  geometry 

Stability  of  EBT  ring-plasma 

Formation  +  tearing  problems  in  RFM,  RF9P,  spheromak 

All  the  old  unsolved  problems  including  bump  on  tail  (  a  particles) . 

Nelson's  cost-hierarchy  of  models,  repeated  at  the  end  of  the  panel  dis¬ 
cussion  write-up,  was  intended  to  persuade  simulators  to  seek  the  most  appro¬ 
priate  model,  matched  to  the  physics  sought,  using  the  simplest  (least  expen¬ 
sive)  first  and  then,  if  need  be,  moving  on  to  more  elaborate  (more  expensive) 
later.  This  ranking  was  not  intended  as  some  form  of  cost-effectiveness  argu¬ 
ment  for  and  against  fluid  versus  particle  simulations. 

Bruce  langdon,  in  Paper  #11,  made  several  strong  suggestions  about 
doing  MFE  computing  more  effectively  and  more  efficiently.  Optimization  pays 
off  very  quickly  (he  gave  experience  with  ZOHAR  showing  that  optimization 
paid  off  almost  as  it  was  done!)  Optimization  must  be  done  selectively,  so 
as  not  to  hurt  in  making  changes  (e.g.,  ZOHAR  particle  boundary  conditions 
are  handled  in  Fortran).  He  advocated:  use  of  monitoring,  intervening, 
postprocessing,  and  linking  to  other  codes;  relying  heavily  on  interactive 
codes,  sharing  the  (large)  data  set  of  the  simulation  code,  doing  consider¬ 
able  computing  and  driving  rapid  graphic  displays.  Postprocessing  is  very 
helpful  in  evoking  good  physical  understanding,  such  as  spectra,  filtering, 
finding  spatial  and  temporal  correlations,  etc.  Fast  graphic  display,  like 
LLL's  THOS,  is  extremely  valuable.  The  pitch  to  MFE  was  not  only  to  update 

software,  but  also  to  update  hardware. 
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V.  RECOMMENDATIONS 

(from  C.  K.  Birdsall's  notes) 

Recommendations  are  both  explicit  and  implicit  in  the  Panel  Dis¬ 
cussion  and  in  the  above.  TO  ALL: 

(f)  in  computing  be  more  effective,  more  efficient; 

(2)  in  physics,  make  results  better  known. 

TO  VOE: 

(1)  be  more  aware  of  the  day-by-day  interactions  within  the  particle 
simulation-theory-experiment  net,  continually  refining  theory  and 
fluid  codes; 

(2)  be  more  aware  that  most  particle  simulation  workers  have  moved 
toward  particle-fluid  models  to  stay; 

(3)  note  that  particle  simulations  are  now  capable  of  working  on  trans¬ 
port  time  scales; 

(4)  plug,  on  your  end,  for  more  efficient  equipment; 

(5)  include  particle  simulation  people  on  your  computer  time  allocations 
commi ttee; 

(6)  reward  inventions,  like  long-time  averaging,  orbit-averaging,  hybrid 
codes,  MHD-particle  codes,  use  of  global  integrals  of  motion,  use 

of  macrocells,  etc.,  when  these  result  in  real  savings; 

(7)  keep  alive  support  of  the  smaller  university  efforts,  which  in 
turn  aid  much  in  keeping  honest  the  larger  national  lab  efforts; 

(S)  penalize  users  with  larger  budgets  who  are  slow  to  use  new  inven¬ 
tion  to  save  time; 

(9)  do  away  with  "if  you  don't  use  it,  you  lose  it"  budget  policy. 
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